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Preface 


My  original  objective  in  performing  this  thesis  was  to  learn  about 
this  strange  fractional  calculus  by  applying  its  principles  to  unsteady 
aerodynamics.  The  original  goal  was  to  analytically  integrate  the 
press ure-downwash  integral  equation  (a  feat  which  had  never  been 
accomplished  before)  by  sinplifying  the  mathematics  through  the  use  of 
the  fractional  calculus.  As  the  research  progressed,  the  unsteady 
aerodynamic  problem  proved  to  be  more  challenging  than  originally 
anticipated  and  the  analytic  integration  was  not  achieved.  However, 
this  treatment  is  the  first  recorded  attempt  to  model  three-dimensional 
unsteady  aerodynamic  forces  on  wings  with  fractional  calculus.  By 
being  the  first,  there  was  a  great  deal  of  ground  to  cover  and  a  fair 
amount  of  dead  ends  discovered.  However,  1  was  able  to  simplify  the 
unsteady  three-dimensional  aerodynamic  problem  by  the  development  of  the 
equivalent  Theodorsen  function.  Hus  development  will  permit  the 
unsteady  aerodynamic  loads  on  a  finite  airfoil  to  be  written  in  a 
compact  mathematical  form  and  easily  evaluated  and  applied  to  control 
system  design.  The  thesis  is  written  in  such  a  way  to  allow  an 
individual  unfamiliar  with  the  material  to  read  and  understand  the 
concepts  and  to  continue  the  work  if  desired. 

I  wish  to  give  a  special  thanks  to  Lt  Col  Ron  Bagley  for  his 
tolerant  listening  abilities  and  strong  guidance.  I  would  also  like  to 
thank  Captain  Oreg  Warhol  a  for  providing  me  with  an  appreciation  of 
mathematics  and  for  his  encouragement  and  understanding  during  a 
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difficult  personal  crisis.  Finally,  I  would  like  to  express  my 
heartfelt  appreciation  to  my  wife,  Suzie,  and  daughter,  Alex,  for  their 
understanding  and  strong  support  in  helping  me  achieve  one  of  my  goals. 

I  would  like  to  dedicate  this  thesis  to  two  special  people. 

First,  to  my  mother  whose  strength  and  determination  through  a  recent 
divorce  inspired  me  to  continue  under  difficult  circumstances.  And  to 
Alex,  on  the  anniversary  of  her  first  birthday,  whose  acccnpl ishments 
will  someday  outshine  my  own. 
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Abstract 


A  fractional  calculus  model  is  developed  for  the  kernel  function  under 
incompressible  subsonic  flow  conditions  for  rectangular  planform 
airfoils  with  snail  aspect  ratio.  A  model  valid  for  restricted  regions 
of  the  kernel  function  for  compressible  subsonic  flow  conditions  is  also 
developed.  Additionally,  a  method  for  numerical ly  solving  the  pressure- 
downwash  integral  equation  for  rectangular  planform  wings  of  aspect 
ratio  two  through  ten  in  incompressible  flow  is  developed.  An 
equivalent  Theodorsen  function  for  three-dimensional  unsteady  flow  is 
developed,  enabling  the  use  of  the  siitpler  two-dimensional  aeroelastic 
equations  of  motion  to  fully  capture  the  more  ccnplicated  three- 
dimensional  effects.  ^  . 
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MODELLING  AND  ANALYSIS  OF  KERNEL  FUNCTION  AND  DEVELOPMENT  OF  EQUIVALENT 
THEODORSEN  FUNCTION  FOR  3-D  AEROELASTIC  ANALYSIS 


I .  Introduction 


This  thesis  presents  the  development  and  results  of  two 
alternative  approaches  to  the  prediction  of  aerodynamic  loads  produced 
by  time  dependent  motions  of  thin  wings  in  rectilinear  subsonic  flight. 
The  first  approach  taken  is  an  attempt  to  directly  integrate  the 
pressure-downwash  integral  equation.  The  second  approach  taken  is  the 
development  of  the  equivalent  Theodorsen  function  for  three-dimensional 
unsteady  aerodynamics. 

The  first  of  the  two  approaches  was  motivated  by  Bagley  (4:16)  who 
demonstrated  the  ability  of  fractional  calculus  to  model  the  three- 
dimensional  kernel  function  at  the  conditions  given  in  (31:718).  The 
three-dimensional  kernel  function  is  the  transfer  function  relating  the 
airloads  to  the  downwash  (vertical  velocity)  for  a  wing.  This  first 
approach  is  an  attempt  to  model  the  transcendental  nature  of  the  kernel 
function  with  a  mathematically  simpler  function.  The  kernel  function  is 
defined  throughout  the  carpi ex  s-plane,  but  because  of  the  carpi icated 
mathematics,  analysis  is  usually  restricted  to  the  imaginary  axis.  This 
restricts  the  use  of  the  kernel  function  to  stability  analyses  such  as 
flutter  and  has  little  benefit  to  the  control  system  designer.  A  simple 
model  which  captures  fully  the  behavior  of  the  kernel  function 
throughout  the  entire  s-plane  could  be  of  valuable  use  in  active  control 
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algorithms.  The  model  of  the  kernel  function  (the  transfer  function) 
could  be  applied  directly  to  control  system  design  to  increase  the 
capabilities  of  active  control  of  unsteady  aerodynamic  loads. 

The  second  of  the  two  approaches  develops  the  equivalent 
Theodors en  function.  The  equivalent  Theodors en  function  is  constructed 
by  relating  two-  and  three-dimensional  lift  coefficients  in  unsteady, 
incompressible  flow.  The  three-dimensional  lift  coefficients  were 
evaluated  numerically  by  a  method  developed  in  this  thesis  which  is  a 
more  general  form  of  a  method  developed  by  Ueda  and  Dowell  (27:350).  A 
fractional  calculus  model  is  used  to  capture  the  behavior  of  the 
equivalent  Theodorsen  function.  The  modelling  was  motivated  by  Swinney 
(24:5)  who  successfully  modelled  the  two-dimensional  Theodorsen  function 
using  fractional  calculus.  The  relatively  simple  form  of  the  equivalent 
Theodorsen  function  permits  the  inclusion  of  the  three-dimensional 
effects  in  the  two-dimensional  equations  of  motion.  Three-dimensional 
theory  produces  more  accurate  results  than  those  of  two-dimensional 
theory.  Therefore,  the  equivalent  Theodorsen  function  in  fractional 
calculus  form  should  prove  to  increase  the  accuracy  in  the  two- 
dimensional  equations  of  motion  without  greatly  increasing  the  effort 
required  to  generate  a  solution. 

This  thesis  is  divided  into  eight  chapters  and  five  appendices. 

The  first  chapter  is  the  introduction.  Chapter  two  discusses  the 
background  of  the  kernel  function,  presents  a  computational  form  for  the 
kernel  function,  and  describes  the  general  solution  methodology  of  the 
three-dimensional  aeroelastic  problem.  The  next  chapter  presents  a 
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brief  sunnary  of  fractional  calculus  and  discusses  information  necessary 
to  understand  the  application  in  this  work.  References  for  further 
study  will  be  given.  Chapter  four  develops  the  fractional  calculus 
modelling  of  the  modified  kernel  function  for  both  compressible  and 
incompressible  flow  conditions.  An  approximate  solution  methodology  for 
the  pressure-downwash  integral  equation  is  developed  in  Chapter  five. 

The  equivalent  Theodorsen  function  is  derived  in  Chapter  six  along  with 
accompanying  fractional  calculus  models.  Chapter  seven  shows  the 
utility  of  the  equivalent  Theodorsen  function  in  fractional  calculus 
form  by  presenting  compact  time-domain  representations  to  typical 
aeroelastic  responses.  The  thesis  concludes  with  Chapter  eight. 

Six  appendices  have  been  included  to  present  additional 
information  without  detracting  from  the  flow  of  the  thesis.  Appendices 
A  and  B  supply  additional  fractional  calculus  models  of  the  modified 
kernel  function.  The  parameter  variations  of  the  fractional  calculus 
models  for  the  modified  kernel  function  are  shown  in  Appendix  C. 
Fractional  calculus  models  of  the  coefficient  of  lift  for  wings  of 
aspect  ratios  between  two  and  ten  is  contained  in  Appendix  D.  The 
variation  of  the  parameters  of  the  equivalent  Theodorsen  function  will 
be  shown  in  Appendix  E.  Finally,  a  presentation  of  the  original 
development  of  the  pressure-downwash  integral  equation  is  included  in 
Appendix  F. 
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II .  Background  and  Definition  of  Kernel  Function 


Kussner  is  generally  credited  with  the  development  of  the 
pressure-downwash  integral  equation  and  the  three-dimensional  kernel 
function.  A  presentation  of  his  development  in  English  is  given  in 
(16:1-28)  and  a  sunrary  of  the  development  is  shown  in  Appendix  F.  His 
integral  equation  relating  pressure  to  downwash  is  the  solution  to  the 
wave  equation  written  in  cartesian  coordinates.  The  dependent  variable 
is  chosen  to  be  Prandtl's  acceleration  potential  that  is  directly 
related  to  the  perturbation  pressure  field  (16:2).  The  wing  is  treated 
as  a  nearly  plane  impenetrabl e  surface  S'  lying  in  the  x-y  plane  as 
shown  in  Figure  1.  It  should  be  mentioned  that  the  convention  of  the  z- 
axis  positive  downward  was  adopted  subsequent  to  the  original 
development  to  cast  the  equations  into  a  form  which  produced  results 
compatible  with  those  of  analytical  flutter  analyses.  The  x-y-z 
coordinate  system  and  the  surface  S’  are  assumed  to  move  with  uniform 
velocity  in  the  negative  x  direction.  The  solution  is  forced  to  be 
unique  by  satisfying  three  conditions.  First,  disturbances  must  vanish 
far  away  from  the  wing  and  its  wake.  Second,  the  perturbation  pressure 
can  only  be  discontinuous  within  the  region  of  the  surface  S'.  Finally, 
the  perturbation  pressure  must  vanish  along  the  trailing  edge  of  the 
surface  S'  to  satisfy  the  Kutt.a  condition.  Assuming  a  harmonic  downwash 
and  satisfying  the  three  conditions  just  mentioned  produces  the  unique 
solution  written  in  equation  (1).  The  differential  equation  can  be 
rewritten  as  an  integral  equation  relating  the  downwash  w(x,y,t)  at  any 
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Figure  1:  Definition  of  Coordinate  System 


point  (x,y)  on  the  wing  to  the  perturbation  pressure  coefficient 
AP(£/tt»t)  at  a  point  ,t| )  on  the  wing. 

w(x,y,  t)  _  f +1 
V  Bit  J -i 


t)  K[M,  Jc,x-5,  s(y-Ti)  ]  d^cfn 


The  variables  in  equation  (1)  are  defined  as  follows 

t  semispan  of  wing 

$je,  $te  coordinates  of  leading  and  trailing 

edges 

M  Mach  number 

V  velocity 

K[M,k,x-(  ,s(y-t| )  ]  kernel  function 


dimensionless  chordwise  variables 
referred  to  root  senrichord,  bh 

dimensionless  spanwise  variables 
referred  to  semispan,  ; 

s 

ratio  of  semispan  to  root  semichord 

*1 

root  semi chord 

Equation  (1)  is  referred  to  as  the  pressure-dcrwnwash  integral  equation 
in  three-dimensional  unsteady  aerodynamics.  The  perturbation  pressure 
coefficient  AP  in  equation  (1)  is  defined  as 


AP(£,ti,  t) 


(2) 


The  downwash  and  pressure  coefficient  distributions  in  equation  (1)  are 
both  assumed  positive  downward.  When  the  analysis  is  restricted  to 
wings  of  rectangular  planform,  then  and  j  are  replaced  with  -1 
and  +1  respectively. 

The  pressure- integral  equation  (1)  is  used  in  the  direct  sense  to 
solve  for  the  unknown  pressure  coefficient  distribution  for  a  given 
planform,  a  known  or  assumed  mode  of  oscillation,  and  a  prescribed  set 
of  stream  conditions.  A  pressure  coefficient  distribution  must  be 
determined  which  satisfies  the  edge  conditions  appropriate  to  the 
planform  and  flow  regime  under  consideration  and  which,  when  multiplied 
by  the  kernel  function  and  integrated  over  the  planform  yields  the 
downwash  distribution  corresponding  to  the  mode  of  oscillation  (32:5). 
The  pressure-integral  equation  (1)  is  used  in  the  direct  sense  to  solve 
for  the  downwash  distribution  for  a  given  planform,  an  assumed  or  given 
pressure  coefficient  distribution,  and  given  stream  conditions. 

Numerical  procedures  are  usually  employed  to  solve  the  pressure- downwash 
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integral  equation  as  analytic  solutions  have  not  yet  been  determined  for 
most  realistic  planforms. 

It  should  be  noted  that  the  integration  of  the  spanwise  variable 
will  require  the  use  of  the  finite  part  of  an  infinite  integral  in  the 
Mangier  sense  (32:13)  to  handle  the  double-pole  singularity  in  y-^  which 
is  evident  in  equation  (5).  The  generalization  to  transient  motion  is 
accanplished  by  forming  the  Laplace  transformed  version  of  equation  (1). 
This  is  accomplished  by  replacing  ik  with  complex  s  (2:5) 

ft(x,y;g  )  _  V  [** 

V  8it  J-i 

(3) 

rlem<"  &P{£,  t);  s  )  k[M,  s,x-C,  s(y-Ti)  ]  dldr\ 

Here  the  circumflex  indicates  a  Laplace  transform  and  s=sb^/V  is  the 
dimensionless  Laplace  parameter. 

The  majority  of  references  addressing  the  kernel  function  have 
been  primarily  concerned  with  numerical  computation  of  the  kernel 
function  and  numerical  integration  of  the  pressure-downwash  integral 
equation.  Analytic,  closed  form  solutions  of  the  pressure-downwash 
equation  (1)  have  been  found  only  a  small  number  of  special  planforms 
and  flow  conditions  (32:2).  One  of  the  first  and  most  often  cited  works 
treating  the  kernel  function  is  (31:703-718).  This  pioneering  work  cast 
the  kernel  function  into  a  form  more  amenable  to  numerical  computation 
and  provided  explicit  relations  for  handling  limiting  cases  such  as 
incompressible,  sonic,  and  steady  flow  conditions.  The  paper  also 
addressed  the  nature  of  the  singularities  present  and  provided  series 
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representations  of  the  kernel  function  within  regions  of  its  domain. 

Two  of  these  authors  later  added  another  author  and  produced  (32:1-21) 
which  presented  a  systematic  procedure  for  solving  the  integral 
equation.  This  method  involved  assuming  a  plausible  pressure  mode  shape 
which  added  more  singularities  and  difficulties  to  the  problem  but  the 
method  does  produces  satisfactory  results.  Next  a  series  of  at  least 
twenty  different  attempts  to  solve  the  integral  equation  reached  the 
literature.  One  of  the  more  popular  is  that  of  Reissner  who  developed  a 
method  based  on  a  lifting  surface  approach  to  the  solution  of  the 
integral  equation  (21:1-39,22:1-97)  that  was  exact  in  the  limit  as  the 
span  approached  infinity  and  resembled  lifting  line  theory  at  steady 
flow.  The  next  series  of  attempts  to  solve  the  integral  equation 
involved  approximate  solutions  based  upon  discrete  element  methods  and 
was  coincident  with  the  development  of  faster  computers.  Among  these 
methods  are  the  vortex  lattice  (20:1-492),  doublet  lattice  (1:279-285), 
doublet  point  (27:348-355),  finite  element  techniques  (12:626-633)  and 
variational  techniques  (29:492-498)  just  to  name  a  few.  Landahl  and 
Stark  (13:2049-2060)  wrote  a  complete  synopsis  of  the  progress  made  in 
the  solution  of  this  problem.  Each  of  these  methods  produces 
satisfactory  results,  some  more  easily  than  others.  Coincident  with 
this  development  was  the  development  for  solutions  to  the  larger  class 
of  problems,  namely  the  non-planar  wing  (14:1045-1046).  The  number  of 
different  methods  and  procedures  attempted  to  obtain  the  best  solution 
are  a  testament  to  the  difficulty  of  obtaining  solutions  to  the 
pressure-downwash  integral  equation. 

I 
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Computational  Form  of  Kernel  Function 

Accurate  computation  of  the  kernel  function  is  essential  if  the 
kernel  function  is  to  be  modeled  with  any  degree  of  accuracy.  The 
kernel  function  poses  some  numerical  impediments  that  mast  be  carefully 
avoided  for  accurate  results.  The  numerical  difficulties  arise 
primarily  from  the  presence  of  both  discontinuities  and  singularities. 
The  kernel  function  can  be  written  as  (32:7) 


K[M,k,xQ,  s(y- r| ) ] 


K*\M,k,x0,  s{y-r\)' 
b%s2(y-x\)3 


(4) 


where  K*  is  denoted  the  modified  kernel  function  and  the  shorthand 
notation  xjj=x-£  is  adopted  for  brevity.  This  form  of  the  modified 
kernel  function  is  particularly  helpful  because  it  shows  the  second 
order  singularity  at  y=r)  that  exists  in  the  kernel  function  and  admits  a 
less  difficult  calculation  of  the  modified  kernel  function,  K*.  The 
modified  kernel  function  is  defined  in  equation  (5)  (32:7)  where  and 

Ii  are  modified  Bessel  functions  of  the  first  order  and  first  and  second 

▲ 

kinds  respectively,  1^  is  the  modified  Struve  function  of  the  first 
order  (8:374-379,495-502)  and  f3=(l-M^)!^.  The  modified  kernel  function, 
K*  does  not  contain  any  singularities  but  it  does  contain  one 
discontinuity  at  Xg=x-{=0  and  yg=y-n=0  which  can  be  shown  to  have  the 
limiting  form  shown  in  equation  (6)  (32:7). 
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K*  [  M,  k,x0,  s(y-rj)  ]  =  Jcs  |y-ti  |  +  Jrs|y-t]|J^(irs|>r-Tj|) 


■jffe  -  Ms/xZ'Pe'iy-V 

+  i-|ics|y-n|[  JiUcsly-T^O -Lidcsjy-t^l)]  +  ------ . - 

2  >/4  +  Vs2(y-T))2 


i  -  fSTT - ;[*«  -  V^  +  p^ty-n)*] 

-  iks\y-r)\J 


/I 


.ei*«|y-T,|tdT 


+  T' 


(5) 


lim  K‘[M,  k,x0,  fl(y-n)] 
y-r|-*0 


2e 


0 


-ik Xq 


*0>0 


XQ<0 


(6) 


This  form  of  the  modified  kernel  function  as  shown  in  equation  (5) 
exists  for  steady  motion  (k=0)  as  well  and  can  be  shown  to  have  the 
following  limiting  form  (32:7,31:710) 


=  1  +  — ■  — -  (7) 

J4  +  sa(y-  n)a 

Although  the  calculation  of  the  modified  kernel  function  may  seem 
trivial,  there  are  two  numerical  problems  hidden  within  the  equation. 
The  first  problem  is  the  quantity  I^-Lj  and  the  second  is  the  definite 
integral  term  in  equation  (5).  The  quantity  I^-Lj  within  the  modified 
kernel  function  is  difficult  to  accurately  compute  for  large  argument 
because  both  individual  functions  grow  unbounded  for  large  arguments. 
This  problem  can  be  eliminated  by  using  the  definition  (33:425) 
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J^JcsIy-nl)  -  ^(Jcsly-Til)  =  l^g \y~~n  1  f1  fil?re-*»\y-i,\* dx  (8) 

1C  J  0 

This  definite  integral  could  be  handled  by  numerical  integration  but 
there  would  be  an  unknown  amount  of  uncertainty  associated  with  the 
result  depending  upon  such  factors  as  step  size  associated  with  the 
integration  routine.  An  alternate  approach  is  to  approximate  the 
integrand  with  a  function  which  is  integrable  in  closed  form  and 
evaluate  the  result.  This  approximation  approach  is  used  in  this  case 
because  the  error  of  the  approximation  is  small.  In  some  instances, 
accurate  approximations  to  the  integrand  may  not  be  available  and 
numerical  integration  would  be  the  only  recourse.  The  result  is  the 
following  expression  (32:8) 

I^iesly-Til)  -  Lj(ics|y-r||)  - 

(9) 


a2  +  a. 


3k2s3  (y-T]} 2 


where 


apl.0085 
8^=1.3410 
a, =1.0050 
aL=0.8675 
85=0.4648 
a5=0.9159 
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This  approximation  introduces  an  error  of  approximately  0.4  percent  in 
the  vicinity  of  ksjy-^|=  4  and  is  less  throughout  the  domain  of  the 
function  (32:8) . 


The  second  numerical  difficulty  is  in  the  evaluation  of  the 
integral  term  of  the  modified  kernel  function,  equation  (5).  This 
integral  has  not  been  solved  in  closed  form  to  date.  Numerical 
integration  once  again  is  possible;  instead,  an  approximation  to  the 
integrand,  integrable  in  closed  form,  exists  (32:8)  and  will  be  used  to 
minimize  the  uncertainty  in  the  result.  The  following  approximation  can 
be  used  to  replace  the  integrand  in  the  integral  of  the  modified  kernel 
function 


- - -  -  l  -  -  a2e"r,T-  a3©"rjTsinwt  (10) 

Vl  + 


where 

ai  =0.101 

32=0.899 
a«=0. 09480933 
r,=0.329 
r2=1.4069 
r3=2.90 


The  maximum  error  of  this  approximation  is  about  0.24%  near  x=1.5  and 
this  function  possesses  the  same  limiting  value  as  the  integrand  at  the 
two  limits  of  integration  (32:8).  However,  this  expression  is  only 
valid  for  positive  values  of  x  so  the  integral  must  be  broken  up  into 
two  regions  and  a  change  of  variables  made  on  one  of  the  two  parts  to 
write  the  complete  integral  in  terms  of  positive  limits  of  integration 
where  the  limits  of  integration  are  defined  as  follows. 
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t 


x0>Ms\y-r\ 


(11) 


L 


a*-d 


vT 


ei*«|y-n)T 


+  X' 


fd  ** - 5 - e‘iks)y_,,ltdx 

Jo  yTTir 


x0<tfs|y-ri|  (12) 


a* 


*o 

pas|  y-1) 


(13) 


d  =  MyJx%  +  $2  a2  (y-i))2  (14) 

pas|y-Ti| 

The  modified  kernel  function  is  now  in  a  form  which  can  be  readily  and 
accurately  evaluated  throughout  the  domain  of  any  wing  planform. 

The  solution  of  the  integral  equation  (1)  involves  the  integration 
of  the  kernel  function.  Although  the  solution  i3  not  specifically 
addressed  in  this  section,  it  is  appropriate  to  discuss  the  special 
integration  required  in  the  sparrise  integration.  The  chordwise  £ 
integration  does  contain  a  finite  jtmp  discontinuity  which  can  be  easily 
handled  by  separating  the  integral  into  two  regions,  solving  each 
individually  by  appropriate  means.  The  s panwise  integration  contains  a 
double-pole  singularity  at  y-tl=0  which  necessitates  the  use  of  the 
Mangier  finite  part  of  an  infinite  integral  (32:13)  as  shown  in  equation 
(15).  The  integration  of  a  function  F(^)  other  than  F(ij)=l  is  more 
difficult  in  general  and  can  only  be  performed  provided  the  function  F 
is  not  singular  at  tj=0.  In  most  instances ,  the  integration  is 
sufficiently  difficult  to  warrant  analysis  by  numerical  means.  Van 
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(15) 


f y<  dt\  =  lim  ry*  +  r**c  dn  _  2 
y-(  (y- T])a  e~°  Jy- c  </y+«  (y-T))2  e 


liml 

\  1 

y-*  +  1  y*t 

ll 

e-0  | 

1  y- nJ 

y-c  y-V  Jy+t 

« 

Nierkerk  (30:1196)  and  Ueda  and  Dowell  (27:350)  have  both  developed 
approximate  numerical  methods  capable  of  accurately  evaluating  such  an 
integral.  The  first  method  by  Van  Nierkerk  (30:1196)  uses  a  special 
Gauss  quadrature  rule  which  is  generally  accurate  with  as  few  as  four 
abscissa  points  for  most  functions.  The  second  method  by  Ueda  and 
Dowell  (27:350)  develops  an  equivalent  expression  for  the  integrand 
under  the  singular  integral  which  can  be  used  in  discrete  element 
solutions  to  the  pressure-downwash  integral  equation. 

General  Solution  Methodology 

This  section  outlines  the  modification  of  the  general  aerodynamic 
problem  into  a  form  suitable  for  solution  by  the  approximate  methods 
described  earlier.  The  pressure-downwash  integral  equation  is  used  in 
the  direct  sense  of  lifting  surface  theory  to  solve  for  the  unknown 
pressure  coefficient  distribution  given  a  known  or  assumed  dewnwash 
distribution.  As  there  are  only  a  very  limited  nunber  of  special  cases 
for  which  closed  form  solutions  are  exist,  the  general  problem  must  be 
solved  by  approximate  nunerical  procedures  as  discussed  in  the  previous 
section.  The  wing  is  assumed  to  be  undergoing  a  displacement  H(x,y,t) 
which  is  represented  as  the  superposition  of  rigid  body  and  elastic 
modes  of  vibration 
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+ 


(16) 


H{x,y,  t)  -  hx  {x, y)  qx  ( t)  +  {x,y)  q2  ( t) 

+  hJ{x,y)qj(t)  +  ... 

Linear  superposition  is  possible  because  the  boundary  conditions  used  by 
Kussner  were  linear,  preserving  the  linearity  of  the  wave  equation 
operator  and  hence  the  solution,  the  integral  equation  (1),  is  itself  a 
linear  operator.  For  sinusoidal  oscillations  we  have  the  following 

qjit)  -  (17) 

where  the  tilde  indicates  the  magnitude  of  the  displacement  in  the  j-th 
mode,  M  is  the  frequency  of  oscillation,  hj(x,y)  describes  the  shape  of 
the  j-th  mode  and  i=(-l)^2.  The  downwash  w(x,y,t)  for  a  given 
displacement  H(x,y,t)  is  given  by  the  following  (32:5) 

wlx.y.t)  U8) 

This  can  be  rewritten  using  the  above  relations  in  equations  (16)  and 
(17)  for  the  j-th  mode  of  vibration. 

it)h,  U.y)Sd£L  (19) 

Seeking  the  downwash  as  the  product  of  a  shape  function  and  an  unknown 
time  dependent  magnitude  reduces  the  expression  above  to  the  following 

^ir1  ■(£  ♦")*»<*•!'»  <20) 

where  the  downwash  has  been  assumed  as 
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(21) 


w(x,y,  t)  =&(x,y)ei*t 
The  pressure  coefficient  distribution  on  the  wing  can  also  be  put  into  a 
similar  form 

AP(5,q,t)  =  Ai5^,^;  ic)  eimt 

The  shape  function  in  equation  (22)  is  shown  with  k,  the  reduced 
frequency,  as  parameter.  By  choosing  the  downwash  mode  shape,  the 
pressure  coefficient  is  forced  to  behave  as  specified  by  the  kernel 
function.  The  pressure  coefficient  is  a  function  of  reduced  frequency 
to  counter  the  frequency-dependence  in  the  kernel  function.  The  type  of 
notation  in  equation  (22)  allows  for  easy  transition  to  the  Laplace 
domain  formulation  in  equation  (3).  The  use  of  these  equations  produces 
the  following  relationship  between  the  shape  of  the  downwash  and  the 
shape  of  the  pressure  distribution  throughout  the  spatial  domain  of  the 
airfoil 


tHx,  y)  _  V  f*1 
V  8kJ-i 

(23) 

;  k)  K  [M,k ,X-Z,  s  (y-Ti )  ]  d£c#n 

For  a  general  problem,  the  general  displacement  shape  is  decomposed  into 
a  finite  rnirber  of  modes  (either  rigid  body  or  elastic  modes).  Each 
mode  shape  is  handled  separately,  solving  the  integral  equation  by  an 
appropriate  means  for  the  pressure  coefficient  distribution  responsible 
for  the  assumed  downwash  distribution.  The  summation  of  the  individual 


modal  pressure  coefficient  distributions  produce  the  total  resultant 


pressure  coefficient  distribution  for  the  general  condition.  The 
pressure  coefficient  distribution  can  in  turn  be  used  to  compute 
generalized  aerodynamic  forces  acting  on  the  wing. 
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Ill ,  Fractional  Calculus 


The  mathematics  of  fractional  calculus  is  nearly  as  old  as  that  of 
classical  calculus.  The  fractional  calculus  treats  derivatives  and 
integrals  of  fractional  order  and  is  not  solely  restricted  to  integer 
order  as  is  the  classical  calculus  (23:115).  An  interesting  discussion 
of  the  historical  development  of  fractional  calculus  has  been  written  by 
Ross  (23:115-122).  The  early  engineering  applications  of  fractional 
calculus  were  precipitated  by  the  observation  that  the  stress  relaxation 
phenomenon  of  viscoelastic  materials  appeared  to  proportional  to  time 
raised  to  the  fractional  power.  This  observation  in  turn  suggested  a 
fractional  order  time  derivative  rather  than  an  integer  order  time 
derivative  in  the  equations  of  motion  for  a  system  comprised  of  the 
material.  This  discovery  spawned  a  renewed  interest  in  the  fractional 
calculus  in  the  twentieth  century.  The  fractional  derivative  can  be 
defined  as  the  inverse  operation  of  fractional  integration  attributed  to 
Reimann  and  Liouville  (3:203) 


1 

ra-«) 


dl*m 
d  t1+- 


x(x) 

Jo  (t-T)m 


dx  ; 


(24) 


0 < a <  1  ,  meN 

One  especially  convenient  feature  of  the  fractional  derivative  is 


lx(fc)]}  =  s—S£(x(t))  <25> 

which  shows  differentiation  in  the  Laplace  domain  is  equivalent  to 
multiplication  by  the  quantity  sa+*  (3:203),  3  being  the  general  Laplace 
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variable.  A  similar  relation  holds  in  the  Fourier  domain  when  s  is 
replaced  by  i*  and  Fourier  transforms  are  used  instead  of  Laplace 
transforms  in  equation  (25).  For  example,  the  Laplace  domain 
representation  of  the  fractional  derivative  of  order  one  half  of  the 
function  f(t)=e'it  is  s^/(s+a).  These  mathematical  features,  coupled 
with  the  apparent  correct  description  of  viscoelastic  phenomenon,  have 
generated  the  renewed  interest  in  fractional  calculus  and  its 
applications  to  engineering  problems.  The  fractional  derivative 
operator  Ef^C  ]  is  a  linear  operator  so  all  the  mathematical 
conveniences  associated  with  linearized  problem®  can  still  be  utilized 
when  the  fractional  calculus  is  included  in  the  problem. 
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IV.  Model  of  the  Modified  Kernel  Function 

The  primary  emphasis  of  this  portion  of  research  is  to  develop  a 
model  which  accurately  describes  the  frequency  dependent  and  chordwise 
and  s panwise  spatially  dependent  behavior  of  the  kernel  function.  It 
should  be  emphasized  that  the  modelling  will  be  performed  on  computed 
values  of  the  kernel  function  in  an  attempt  to  simplify  the  mathematics 
involved  in  its  analysis  and  is  not  performed  to  actual  test  data.  One 
major  goal  is  to  substitute  relatively  simple  functions  for 
transcendental  ones,  reducing  the  effort  presently  required  to  calculate 
the  kernel  function  along  the  imaginary  axis  and  throughout  the  complex 
s-plane.  Additionally,  if  a  suitable  model  could  be  found  which  was 
more  amenable  to  analytic  integration  than  the  kernel  function  itself, 
equations  of  motion  for  the  finite  airfoil  could  be  written  directly  as 
functions  of  this  integration  and  evaluated  directly  in  closed  form. 
Finally,  given  the  first  two  goals,  if  the  model  is  of  a  form  which  has 
a  closed  form  inverse  Laplace  transform,  then  this  transform  can  be  used 
as  the  transfer  function  in  control  system  analyses.  The  kernel 
function  can  be  thought  of  as  a  transfer  function  spatially  convolved 
with  the  pressure  coefficient  distribution  to  produce  the  downwash 
distribution.  A  model  which  describes  the  behavior  of  the  kernel 
function  along  the  imaginary  axis,  integrable  in  closed  form,  and 
available  for  analytic  continuation  into  the  complex  plane  is  the  target 
of  the  project.  Given  this  descriptive  model  can  be  developed,  a 
tractable,  analytically  integrable  transfer  function  that  can  be 


spatially  convolved  to  produce  the  response  of  the  system  to  general 
input  will  be  obtained. 

Model  Development 

The  candidate  model  must  be  capable  of  capturing  the  variation  of 
the  kernel  function  over  a  wide  range  of  reduced  frequencies  along  the 
inaginary  axis  with  the  ability  to  be  analytically  continued  into  the 
entire  complex  plane  without  difficulty.  The  range  of  reduced  frequency 
has  to  be  sufficiently  large  to  insure  that  both  the  small  and  large 
argument  asymptotic  behavior  of  the  kernel  function  was  captured.  As 
shown  earlier,  the  kernel  function  has  the  second  order  singularity 
which  could  cause  problems  in  developing  an  accurate  model,  especially 
in  the  vicinity  of  the  singularity.  Therefore,  the  modified  kernel 
function  shown  in  equation  (5)  earlier  will  be  modelled  as  it  is  more 
well  behaved  mathematically.  It  should  be  noted,  however,  that  the 
singular  nature  of  the  kernel  function  shown  in  equation  (4)  can  be 
included  before  the  integration  is  performed. 

The  modified  kernel  function  will  be  modelled  with  the  functional 
form  shown  in  eqaution  (26).  Kg*  is  given  by  the  steady  value  computed 
using  equation  (7).  This  equation  will  be  shown  later  in  Chapter  VII  to 
be  comprised  of  two  functions  which  are  derivatives  of  two  general  order 
Mittag-Leffler  functions  (17:102)  with  step  functions  as  the  leading 
coefficients.  Bagley  (5:742)  has  shown  a  model  similar  in  form  is 
appropriate  for  describing  the  frequency  dependent  behavior  of  the 
modulus  of  a  viscoelastic  material.  In  Bagley's  model  (5:742),  a 
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K'  [M,k,x0,y0] 


Qx{M,x0,y0)  (ijc)e»(*>Wa) 
1  +Qi(M,x0ly0) 


05  (M,x0,y0) 

1  +07  (M,x0,y0)  (ik)  0*<A''x°'y‘,) 


(26) 


single  Mittag-Leff ler  function  was  sufficient  to  capture  the  behavior  of 
the  modulus.  Swinney  (25:5)  was  able  to  capture  fully  the  frequency 
dependent  behavior  of  the  Theodorsen  function  for  two-dimensional  flow 
(25  :418)  using  a  single  Mittag-Leff ler  function.  The  Theodorsen 
function  relates  the  variation  of  circulatory  lift  to  downwash  of  a  flat 
plate  undergoing  simple  harmonic  motion  (25:1).  Swinney's  model  was 
also  extended  to  include  the  laplace  variable  for  arbitrary  motion.  The 
success  in  modeling  the  two-dimensional  aerodynamic  function  led  Bag  ley 
(4:16)  to  suppose  the  three-dimensional  kernel  function  could  be 
modelled  with  a  similar  functional  form.  Early  investigation  over  a 
broad  band  of  reduced  frequencies  demonstrated  Bagley's  form  of  model 
(4:16)  was  incapable  of  capturing  the  behavior  in  the  higher  subsonic 
region  of  the  kernel  function  (M=0.8+) .  The  model  shown  in  equation 
(26)  with  two  terms  rather  than  one  was  adopted  in  an  attempt  to  capture 
better  the  properties  of  the  modified  kernel  function  over  a  wide  range 
of  reduced  frequencies.  The  fractional  calculus  based  model  is 
especially  convenient  because  the  analytic  continuation  is  automatic 
with  the  substitution  of  the  general  Laplace  parameter  for  the  purely 
imaginary  argunent.  Additionally,  this  type  of  functional  relationship 
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is  capable  of  capturing  two  independent  phase  lags  because  of  the  two 
independent  derivatives  of  the  Mittag-Leffler  functions  present. 

There  is  actually  seme  empirical  foundation  for  a  model  of  this 
form.  In  two-dimensional  flow,  the  wing  is  considered  infinite  in  span 
and  the  only  wake  produced  is  along  the  trailing  edge.  Under  steady 
flow  conditions,  the  wake  leaves  the  trailing  edge  smoothly.  As  the 
reduced  frequency  is  increased,  the  wake  no  longer  leaves  smoothly  but 
instead  produces  disturbances  along  the  trailing  edge  which  propagate 
forward  and  change  the  effective  angle  of  attack  of  the  wing.  For  wings 
of  finite  span,  the  wing  tips  may  also  create  disturbances  independent 
of  those  along  the  trailing  edge  which  cause  additional  disruption  in 
the  motion  of  the  airfoil.  This  type  of  phenomenon  may  require  two 
functions  with  i„.  separate  phase  lags.  The  functions  in  equation  (26) 
are  capable  of  representing  two  separate  phase  lags  because  the  8j  and 
8j  valuer  are  independent,  one  capturing  the  trailing  edge  effects  and 
the  other  capturing  the  wing  tip  effects. 

The  model  shown  in  equation  (26)  is  a  direct  function  of  reduced 
frequency  with  the  parameters  being  functions  of  Mach  number  and 
dimensionless  chordwise  and  s panwise  spatial  variables .  The  Mach 
numbers  of  interest  ranged  from  inccnpressible  (M=0)  up  to  high  subsonic 
(M=0.8+).  The  range  of  the  dimensionless  chordwise  variable  would  have 
to  include  values  ranging  from  Xgc[-2,2].  The  range  of  the 
dimensionless  spanwise  variable  would  depend  upon  the  aspect  ratio  of 
the  given  airfoil  syflc[-2s,2s]  where  s  is  the  ratio  of  semispan  to  root 
semi chord.  A  general  comprehensive  model  which  could  describe  the 
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entire  range  of  the  modified  kernel  function  was  not  discovered.  The 
large  number  of  independent  variables  created  additional  difficulties  in 
determining  a  suitable  functional  form  for  the  model  parameters,  . 

j 

However,  many  regions  were  adequately  modelled  and  all  regions  were 
fully  investigated  and  analyzed. 

As  evidenced  by  the  form  of  the  model  in  equation  (26)  and  the 
nature  of  the  modified  kernel  function  discussed  in  equation  (5),  this 
curve  fitting  is  a  highly  nonlinear  problem  and  a  suitable  method  of 
nonlinear  optimization  had  to  be  found  to  solve  the  problem  in  a  timely 
manner.  For  models  with  a  small  nunber  of  parameters,  it  is  sometimes 
sufficient  to  perform  manual  regression.  However,  in  this  case  with 
eight  independent  parameters,  a  nonl inear  least  squares  regression 
routine  using  the  modified  Levenberg-Marquardt  algorithm  (18:431-441)  to 
generate  a  sequence  of  approximations  to  the  minimum  point  was  used  to 
expedite  the  process.  This  algorithm  uses  a  "trust  region"  approach 
with  a  bounded  step.  The  Jacobian  is  needed  to  optimize  the  parameters 
of  the  model  and  is  computed  numerically  in  this  instance  by  forward 
finite  differences  (10:243).  The  IMSL  routine  RNLIN  which  implements 
this  algorithm  was  used  (10:239)  in  the  present  analysis.  There  is  a 
certain  amount  of  insight  required  to  use  such  a  routine  in  this  type  of 
application. 

This  routine  is  designed  for  regression  of  real  valued  functions, 
not  complex  valued  functions.  In  performing  a  modelling  of  the  modified 
kernel  function  a  model  must  be  constructed  which  adequately  predicts 
the  behavior  of  both  the  real  and  imaginary  parts  of  the  function.  A 
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simultaneous  regression  performed  on  both  the  real  and  imaginary  parts 
is  desired.  Unfortunately,  in  this  present  work,  the  simultaneous 
regression  was  not  readily  programmable  and  was  not  used. 

Two  different  attempts  were  made  to  trick  the  routine  into 
performing  a  simultaneous  regression.  The  routine  requires  a  real¬ 
valued  function  which  computes  the  error  between  the  function  and  the 
model  at  a  discrete  number  of  points.  This  is  the  only  possible  input 
to  the  routine  which  could  potentially  be  modified  to  achieve  a 
simultaneous  fit.  The  first  attempt  was  to  return  to  the  routine  the 
magnitude  of  the  complex  error  between  the  points.  The  real  part  caused 
the  routine  to  diverge  in  this  instance.  Another  attempt  was  made  to 
treat  the  real  and  imaginary  parts  as  one  real -valued  curve  by 
translating  the  imaginary  part  to  begin  at  the  tail  of  the  real  part. 
This  failed  as  well  because  the  derivative  of  the  curve  was,  in  general, 
discontinuous  at  the  point  where  the  two  parts  were  joined.  Without  a 
quick  method  of  performing  simultaneous  regression,  the  regression  was 
performed  solely  on  the  imaginary  part.  The  imaginary  part  proved  to  be 
more  amenable  to  fitting  than  the  real  part.  When  the  real  part  of  the 
model  (based  upon  the  regression  of  the  imaginary  part)  was  compared  to 
the  real  part  of  the  modified  kernel  function  it  was  generally  observed 
to  produce  a  satisfactory  fit  as  well.  Normally,  the  real  part  of  a 
complex  function  would  not  be  well  modelled  by  specifying  the  imaginary 
part.  However,  given  a  linear  transfer  function,  the  imaginary  and  real 
parts  can  be  shown  to  be  related  through  derivative  operators  using 
Fourier  transforms.  In  this  work,  a  good  model  for  the  imaginary 
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produced  a  good  model  for  the  real  part  as  well,  evidence  of  the  strong 
linear  nature  of  the  modified  kernel  function. 

Another  obstacle  in  this  IMSL  routine  is  the  requirement  for  a 
starting  point  in  the  optimization  space.  A  good  initial  point  would 
produce  good  regression  results  and  a  bad  initial  vector  would  produce 
no  regression  result  and  unfortunately  no  information  regarding  how  to 
improve  the  initial  value.  The  real  difficulty  in  the  procedure  is 
determining  the  correct  initial  vector;  because  this  function  has  not 
been  modelled  in  this  manner  before,  the  selection  of  the  initial  trial 
vector  was  something  gained  by  experience  and  practice.  Cfcie 
recommendation  for  further  investigation  in  this  area  is  the  use  of  a 
complex  regression  algorithm  for  evaluating  models  of  the  modified 
kernel  function. 

There  are  several  advantages  associated  with  this  IMSL  routine  as 
well.  First,  as  the  routine  is  available  internationally  and  used 
extensively,  the  routine  probably  has  been  extensively  tested  against 
established  test  cases  and  has  no  errors.  Creating  a  routine  takes  a 
large  amount  of  time  to  code  and  test.  Another  advantage  is  speed.  The 
IMSL  routines  are  coded  for  efficient  operation.  The  routine  was  one  of 
the  better  nonlinear  regression  routines  available  in  either  pre-coded 
or  algorithm  format. 

One  interesting  and  helpful  result  discovered  during  the  modelling 
process  was  the  nature  of  the  regression  parameters,  8j .  Given  a  good 
regression  for  one  particular  Mach  nunber,  x^  and  yfl,  a  reasonable 
initial  value  for  another  Xg  and  y5  value  for  the  same  Mach  nimber  could 
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be  interpolated  quite  readily  .  This  resulted  in  family  of  regression 
curves  for  the  various  parameters  of  the  model.  Unfortunately,  this 
behavior  was  not  observed  fc*.  all  regions  of  the  modified  kernel 
function  using  this  form  of  model,  possibly  suggesting  the  present  form 
of  model  is  not  quite  appropriate  in  the  region.  The  existence  of  the 
strong  relationship  between  the  parameters  of  the  model  in  many  regions 
suggests  the  existence  of  seme  type  of  fractional  derivative  model  that 
could  adequately  represent  the  complicated  and  highly  intractable 
modified  kernel  function. 

The  goodness  of  fits  of  the  fractional  calculus  models  were 
minimized  using  a  least  squares  type  of  error  calculation.  The  error  is 
calculated  using  the  following 

f  p»wD  1 

=  j  Elf'V-Wf  <27) 

"D  (.  P-0  J 

where  f(Xp)  is  the  value  of  the  modified  kernel  function  at  the  p-th 
point,  fc(Xp)  is  the  value  of  the  fractional  calculus  model  at  the  p-th 
point  in  the  interval ,  Ng  indicates  the  total  number  of  points  in  the 
interval ,  and  p  is  an  index  locating  the  point  of  comparison.  The 
errors  for  the  models  will  be  given  when  the  model  is  presented.  The 
number  of  points  used  to  compute  the  error  varies  in  the  thesis.  The 
models  and  functions  throughout  the  thesis  were  compared  at  intervals  of 
5  percent  of  reduced  frequency  (the  independent  variable).  Hence,  for 
an  interval  between  zero  and  one,  there  would  be  twenty- one  points 
compared  and  for  an  interval  between  zero  and  three,  there  would  be 
sixty-one  points  of  comparison  between  the  target  points  and  the  model. 
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Model  Results 


The  results  presented  are  restricted  to  thin  rectangular  planform 
airfoils.  It  should  also  be  pointed  out  that  the  modelling  of  the 
modified  kernel  function  was  performed  to  calculated  values  and  not. 
actual  test  data. 

Ccmpressibl e  Flow.  The  modelling  of  the  modified  kernel  function 
was  generally  unsuccessful  for  ccrrpressible  flow  conditions  using  the 
form  of  model  in  equation  (26).  There  were  a  number  of  obstacles 
encountered  in  attempting  to  model  the  modified  kernel  function.  The 
first  attempt  at  modelling  was  performed  at  Mach  number  0.8.  This 
proved  to  be  extremely  difficult  with  the  form  of  model  shown  in 
equation  (26).  The  model  would  predict  behavior  well  in  regions  of  Xg 
between  0.7  and  1.0;  in  other  regions,  the  overall  trends  of  the 
modified  kernel  function  were  followed  but  the  actual  functional  values 
were  generally  not  achieved.  One  possible  reason  for  the  difficulty  is 
the  compressibility  effects  beginning  to  dominate  at  the  higher  Mach 
numbers.  As  the  Mach  number  approaches  one,  shock  waves  mil  begin  to 
form  on  the  wing.  Shock  waves  are  viewed  as  a  non-linear  phenomenon  and 
as  such,  the  kernel  function,  which  is  a  linear  operator,  would  probably 
have  difficulty  capturing  these  effects.  It  can  be  shown  that  the 
modified  kernel  function  at  Mach  number  1.0  (sonic)  is  equal  to  zero 
(xg<0,  all  syg)  compared  to  the  subsonic  case  where  the  modified  kernel 
function  was  equal  to  zero  for  (xg<0,  syg-  0)  equation  (6).  It  is 
postulated  that  possibly  a  third  function  would  need  to  be  added  to  the 
two  function  model  to  help  capture  the  shock  formation  phenomenon.  The 
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investigation  at  Mach  number  0.8  was  abandoned  for  lack  of  conplex  non¬ 
linear  regression  software  to  perform  the  statistical  regression  on  a 
increasingly  difficult  problem.  It  was  evident  at  this  point  that  the 
global  modified  kernel  function  model  would  not  be  obtained.  Another 
Mach  minber  in  the  compressible  regime  was  analyzed  in  an  attempt  to 
characterize  the  modified  kernel  function  at  least  one  Mach  number. 

The  Mach  number  was  reduced  to  0.5  in  an  attempt  to  decrease  the 
dominance  of  the  compressibility  effects,  but  still  retain  enough 
compressibility  to  attempt  to  capture  the  subdominant  effects.  The 
modelling  at  this  Mach  number  produced  same  very  good  results.  The 
modified  kernel  function  was  best  modelled  in  the  regions  with  small  Xg, 
small  syg  and  large  Xg,  large  syg.  Acceptable  models  were  obtained  in 
the  vicinity  of  small  Xg,  large  syg  and  poorer  fits  were  found  for  large 
Xg,  small  syg.  Negative  Xg  values  were  not  analyzed  for  this  case 
because  the  modelling  was  determined  to  be  insufficient  to  capture  fully 
all  of  the  behavior. 

The  resulting  model  for  two  of  the  regions  will  be  shown  in  the 
main  body  of  the  text  and  samples  from  the  remaining  two  regions 
contained  in  Appendix  A.  A  sample  comparison  between  the  real  and 
imaginary  parts  of  the  model  and  modified  kernel  function  in  the  small 
Xg,  small  sy0  region  is  shown  in  Figure  2.  The  model  is  shown  below 
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where 


Kj =1.917 

0,  =-0.027  02=1.205 

03=0.154  0,  =0.094 

©5  =0 . 2  3  6  06=1.O1 

@3=0.064  0j=O.  95 

The  error  for  the  real  and  imaginary  parts  is  2.6xlO"3  and  1.3xl0‘3 
respectively.  The  steady  value  is  computed  using  the  relationship  in 
equation  (7).  As  shown  in  the  figure,  there  is  excellent  correspondence 
between  the  model  and  the  modified  kernel  function  for  both  real  and 
imaginary  parts.  In  this  region,  the  exponential  limiting  form  of 
equation  (6)  is  dominant  and  this  behavior  is  evidenced  in  the  figures 
shown.  If  the  curves  were  extended  to  include  values  of  negative 
frequency,  the  beginning  of  a  finite  jump  discontinuity  could  be  seen, 
resembling  a  travelling  wave  at  this  point.  As  mentioned  earlier,  other 
values  of  the  modified  kernel  function  in  the  vicinity  of  this  chordwise 
and  spanwise  location  can  be  interpolated  with  a  certain  amount  of 
accuracy.  There  appear  to  be  a  substantial  number  of  local  minima  in 
the  nonlinear  least  squares  minimization  function  created  by  the  model 
and  the  modified  kernel  function.  Therefore,  the  modelling  of  one 
region  of  the  function  does  not  necessarily  produce  good  results  for 
another  region.  The  compressibility  effects  are  thought  to  be 
contributing  to  the  large  number  of  local  minima  present.  A  true 
complex  regression  algorithm  might  be  more  beneficial . 

Another  region  investigated  was  small  Xj  and  large  syQ .  The  model 
had  increasing  difficulty  capturing  the  behavior  of  the  modified  kernel 
function  for  larger  syg  primarily  due  to  the  more  oscillatory  nature  of 
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REDUCED  FREQUENCY 


Figure  2:  Model  of  K?[0.5,0.1,0.05] 


the  kernel  function  for  large  arguments  (26:174).  A  sample  of  the  model 
in  the  larger  syj  domain  is  shown  in  Figure  3.  The  model  has  the  form 
shown  in  equation  (29)  with  the  real  and  imaginary  parts  having  errors 
of  3. 2x10’ 2  and  7.1xlO'3  respectively.  As  the  reduced  frequency  is 
increased,  the  modified  kernel  function  becomes  oscillatory  which  would 
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REDUCED  FREQUENCY 


Figure  3:  Model  of  K?[0. 5,0. 1,1.0] 
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require  successively  more  terms  in  the  model.  This  may  suggest  an 
infinite  series  with  an  argument  consisting  of  a  function  similar  to  the 
one  used  in  the  model.  This  is  an  area  that  could  be  further  developed. 
The  large  k  behavior  is  important  because  the  large  asymptotes  must  be 
sufficiently  modelled  to  have  an  accurate  inverse  Laplace  transform  (the 
Laplace  transform  and  its  inverse  are  defined  over  the  infinite  complex 
frequency  plane) . 

The  region  of  large  Xg  and  small  syD  is  shown  in  Figure  17  located 
in  Appendix  A.  Once  again,  in  this  region  the  exponential  limiting  term 
is  dominant  and  the  modified  kernel  function  is  highly  oscillatory, 
causing  difficulty  in  the  modelling.  Although  the  imaginary  part  is 
following  the  function  quite  nicely,  the  real  part  is  losing  its  ability 
to  describe  the  function.  A  true  simultaneous  regression  of  both  real 
and  imaginary  parts  of  the  model  might  have  overcome  this  deficiency. 

The  region  comprising  both  large  Xg  and  large  syg  was  handled 
reasonably  well  for  the  particular  values  investigated  but  the 
performance  of  the  model  will  degrade  as  the  reduced  frequency  is 
increased  because  of  the  oscillatory  nature  of  the  modified  kernel 
function  mentioned  earlier.  Figure  18  in  Appendix  A  show  the  behavior 
in  this  region. 

Overall,  the  modelling  of  the  modified  kernel  function  for  Mach 
number  0.5  was  not  completely  satisfactory  and  did  not  provide  a 
complete  modelling  over  the  domain  of  the  wing.  Further  analysis  into 
the  possibility  of  analytic  integration  of  the  integral  equation  as 
originally  planned  was  not  performed  because  of  time  constraints. 


However,  the  modelling  and  investigation  do  give  promise  that  such  a 
model  could  exist,  given  the  correct  form  could  be  determined.  As 
mentioned  earlier,  a  complex  non-linear  regression  routine  would  greatly 
increase  the  ability  to  explore  other  model  forms.  Lack  of  total 
success  in  the  compressible  domain  led  to  the  investigation  of  the 
modified  kernel  function  for  incompressible  flow.  The  supposition  that 
compressibility  of  any  magnitude  greatly  increased  the  complexity  of  the 
problem  will  be  tested  in  the  next  section. 

Incompressibl e  Flow.  The  incompressible  case  would  require 
another  model.  The  model  given  in  equation  (26),  when  tried  in  this 
case,  had  too  many  parameters  for  the  regression  to  be  performed  with 
any  degree  of  consistency.  The  following  model  was  adopted  for  the 
incompressible  flow  conditions  where,  once  again,  the  Kq*  value  is  given 


K^k.x^.y^]  •  Kf  - 


(30) 


1  +»3(x0,y0)  (ilc) 

by  the  steady  limiting  form  in  equation  (7)  and  the  Mach  number,  M,  has 
been  eliminated  from  the  argument  for  brevity.  A  more  thorough  study 
was  conducted  on  the  incompressible  case  than  on  the  other  compressible 
cases  previously  conducted  because  the  incompressible  case  provided  the 
most  opportunity  for  ccnplete  modelling  of  the  modified  kernel  function. 
One  model  describing  the  modified  kernel  for  incompressible  flow 
conditions  throughout  its  entire  domain  was  not  found  either.  However, 
the  model  in  equation  (30)  shows  a  very  strong  potential  for  modelling 
wings  of  small  aspect  ratios  (aspect  ratios  less  than  one-half)  such  as 
those  found  in  the  late  stages  in  turbomachinery.  Hie  modelling  of  the 


modified  kernel  function  in  the  incompressible  regime  discussed  in  this 
section  was  shown  to  be  very  capable  of  describing  the  snail er  aspect 
ratio  wings  throughout  the  domain  of  the  airfoil.  However,  the 
asymptotic  behavior  of  the  model  with  respect  to  the  modified  kernel 
function  for  snail  and  large  asynptotic  values  of  reduced  frequency  has 
not  established  and  is  currently  under  investigation.  Ueda  (26:169- 
174,28:346-347)  provides  asynptotic  expansions  of  the  kernel  function 
for  both  small  and  large  arguments.  These  behaviors  need  to  be  verified 
before  the  model  can  be  said  to  fully  describe  the  modified  kernel 
function  for  small  aspect  ratio  wings.  Samples  of  the  model  in  specific 
regions  are  presented  here  and  in  Appendix  B. 

A  sample  of  the  model  for  small  syg  and  large  Xg  is  shown  in 
Figure  4.  The  model  has  the  following  representation 

**[*,-0.1, 0.5]  =  <  -  — -  (31) 

1  +  03  (ik)  ®4 

where 

Kg*=l.  9987 

0i=1.729  87=0.81 

03=0.267  04=1.34 

with  errors  of  9.9xlO‘J  and  2.4xlO'J  for  the  real  and  imaginary  parts. 

The  agreement  between  the  model  and  the  computed  modified  kernel 
function  is  much  better  than  the  same  case  for  the  Mach  nunber  0.5  case, 
indirectly  supporting  the  supposition  that  compressibility  effects  cause 
difficulty  in  the  modelling  of  the  modified  kernel  function.  This 
particular  case  represents  the  worst  agreement  between  model  and 
modified  kernel  function  for  the  incompressible  case.  Although  not 
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REDUCED  EREQUENCY 


Figure  4:  Model  of  K*[0, 1.0, 0.05] 

perfect,  it  is  reasonable  considering  the  relative  time  required  to 
evaluate  equation  (30)  and  equation  (5).  The  model  also  worked 
reasonably  well  for  large  Xg  and  large  syfl.  The  small  Xg  and  large  syg 
behavior  was  modelled  very  well  using  equation  (30)  as  shown  in 
Figure  5.  The  model  for  this  particular  case  is  written  in  equation 
(32)  with  errors  of  l.lxlO'2  and  4.7xlO'3  respectively  for  the  real  and 
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REDUCED  FREQUENCY 


Pigure  5:  Model  of  K*[0,0. 1,1.0] 


imaginary  parts. 
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Cornparing  this  figure  to  the  corresponding  case  at  M=0.5  (Figure  3) 
shows  the  higher  Mach  number  case  has  a  much  sharper  curve  than  the 
respective  incompressible  case,  demonstrating  the  more  complex  behavior 
due  to  compressibility  effects.  In  addition  to  investigating  the 
positive  Xj  behavior,  negative  Xg  behavior  was  also  analyzed.  The 
behavior  for  small  values  of  negative  Xg  and  small  syg  is  given  in 
Figure  19  which  is  included  in  Appendix  B.  The  model  clearly  captures 
the  behavior  of  the  modified  kernel  function  in  this  region.  Note  that 
t.hsse  curves  will  be  forced  to  zero  as  the  syA  values  are  decreased 
toward  zero  in  keeping  with  the  limiting  form  given  in  equation  (6).  In 
this  region,  where  the  dominant  behavior  is  characterized  by  a  step 
function  of  magnitude  two,  it  is  important  to  include  some  negative 
reduced  frequency  values  in  the  regression  process  to  insure  that  the 
step  behavior,  the  small  argument  asymptotic  behavior,  is  fully 
captured.  The  advantage  of  the  fractional  derivative  model  is  that  it 
is  capable  of  representing  a  step  function  type  of  behavior  quite 
easily.  This  particular  type  of  behavior  in  the  modified  kernel 
function  has  just  been  initial 1}  investigated  but  is  not  presented  here. 
Figure  20  shows  an  example  of  the  model  in  this  region.  The 
difficulties  associated  with  the  modelling  of  the  modified  kernel 
function  occurred  in  the  regions  of  the  step  function  and  the  region 
containing  negative  Xg  and  large  values  of  syg.  As  stated  earlier,  the 
restriction  to  values  of  syg  less  than  one  merely  suggests  the 
suitability  of  the  model  for  smaller  aspect  ratio  wings.  Another  form 
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of  model  might  be  necessary  to  successfully  model  rectangular  airfoils 
of  higher  aspect  ratio. 

It  is  possible  to  construct  curves  describing  the  dependence  of 
the  parameters  of  the  model  upon  Xq  and  syc.  A  sanple  of  the  variation 
of  the  parameters  of  the  model  for  the  incompressible  case  is  contained 
in  Appendix  C.  These  curves  can  in  turn  be  fit  with  polynomials  of  up 
to  order  five  quite  easily.  However,  rather  than  viewing  these  as 
polynomial  curves,  it  is  more  convenient  to  view  the  terms  of  the 
polynomial  models  as  the  first  few  terms  in  an  infinite  power  series  of 
same  unknown  function  which  would  fully  describe  the  behavior  of  the 
parameters.  If  the  infinite  series  can  be  shown  to  represent  a 
mathematically  convenient  functional  expression,  it  could  facilitate  the 
direct  analytic  integration  of  the  integral  equation  or,  at  least,  allow 
for  an  easier  inverse  Laplace  transformation.  These  ideas  have  been 
initially  studied  in  this  work  and  continued  work  in  this  vein  is 
planned  to  determine  if  this  is  a  viable  approach  to  the  solution  of  the 
problem. 

The  modelling  of  the  modified  kernel  function  with  fractional 
calculus  types  of  functions  has  promise.  It  has  been  demonstrated  that 
models  can  be  developed  which  are  sufficient  for  describing  much  of  the 
behavior  for  moderate  Mach  number  but  not  totally  descriptive  over  all 
the  domain  of  the  kernel  function.  Higher  Mach  nutber  behavior  is  more 
difficult  to  capture  but  some  behavior  was  captured  with  a  very  sinple 
model  suggesting  further  investigation.  The  modified  kernel  function 
for  incompressible  flow  was  found  to  be  modelled  quite  well  for  small 
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aspect  ratio  rectangular  airfoils  by  a  fractional  calculus  model  with  as 
few  as  four  parameters.  Although  the  model  asymptotes  for  small  and 
large  values  of  reduced  frequency  have  not  yet  been  veri  result 

that  the  majority  of  the  domain  of  the  modified  kernel  function  could  be 
captured  with  a  simple  fractional  calculus  makes  an  impo’-tu^t  thrust 
into  a  new  area  of  application  for  fractional  calculus. 

One  additional  topic  which  must  be  addressed  before  the  model  can 
be  deemed  complete  is  the  comparison  with  the  modified  kernel  function 
in  the  carrplex  s-plane.  Ashley  (2:7-8)  provide  forms  of  the  modified 
kernel  function  suitable  for  evaluation  throughout  the  complex  plane, 
one  of  which  has  been  used  here  as  equation  (5).  If  this  model  is  shewn 
to  predict  the  behavior  of  the  modified  kernel  function  in  the  carpi ex 
s-plane,  then  the  model  will  be  suitable  for  control  analyses,  not  just 
stability  analyses. 

This  particular  work  was  only  successful  in  modelling  wings  of 
very  small  aspect  ratios  under  one-half.  The  goal  of  finding  a  directly 
integrable  model  of  the  kernel  function  does  seem  unlikely  but  the 
development  of  a  model  that  can  be  inverted  in  the  Laplace  or  Fourier 
sense  seems  very  feasible.  This  type  of  model  would  be  of  extreme 
importance  as  it  would  permit  the  incorporation  of  the  model  of  the 
kernel  function  into  the  equations  of  motion  of  the  airfoil.  If  the 
model  can  be  extended  analytically  into  the  left  half  s-plane,  then 
equations  of  motion  for  general  motion,  not  just  stability,  could  be 
solved  using  the  simple  model  of  the  kernel  function.  There  is  a  great 
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need  for  more  work  in  this  area  as  this  was  the  first  attenpt  at 
performing  such  a  modelling  an  this  function. 
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V.  Approximate  Solution  of  the  Pressure-Downwash  Integral  Equation 

Although  the  approximate  solution  of  the  integral  equation  (3)  or 
(23)  was  not  a  direct  goal  of  the  thesis,  it  became  necessary  after  the 
lack  of  success  experienced  in  Section  IV  discussed  earlier.  The 
approximate  method  is  used  to  calculate  the  inverse  kernel  function  K‘* 
by  sampling  the  integral  equation  and  reducing  the  problem  to  a  linear 
system  of  equations.  The  inverse  kernel  function  can  than  be  used  to 
calculate  numerically  the  generalized  aerodynamic  forces  which  will  be 
used  to  determine  the  equivalent  Theodorsen  function  for  three- 
dimensional  aerodynamics. 

This  chapter  outlines  the  development  of  a  method  to  evaluate 
numerically  the  integral  equation  (3)  or  (23)  which  will  collapse  down 
to  a  form  similar  to  the  doublet  point  method  developed  by  Ueda  and 
Dowell  (27:348-355).  The  doublet  point  method  (27:348-355)  is  a 
variation  of  the  doublet  lattice  method  (1:279-285)  which  is  easier  to 
apply.  Although  the  doublet  point  method  was  not  used  directly,  the 
result  of  the  author's  method  reduces  to  a  form  which  is  close  to  that 
of  the  doublet  point  formulation  warranting  the  reference  and  credit  for 
the  original  work. 

The  procedure  involves  dividing  the  wing  into  a  finite  nunber  of 
panels  and  representing  the  entire  pressure  coefficient  and  downwash 
distributions  over  the  small  panel  by  a  single  concentrated  load  with  a 
yet  to  be  determined  magnitude.  This  is  mathematically  equivalent  to 
treating  the  pressure  coefficient  distribution  as  a  series  of 
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acceleration  doublet  points  (27:349)  represented  as  Dirac  delta 
functions .  A  down  wash  mode  is  assured  and  a  1  inear  system  of  equations 
is  formed  by  the  sampling  caused  by  the  delta  functions.  The  linear 
system  is  solved  to  determine  the  unknown  coefficients  of  the  delta 
functions  which  define  the  pressure  coefficient  distribution  in  a 
discrete  sense. 

The  integral  equation  (23)  is  separated  as  shown  in  equation 
(33)  to  permit  special  treatment  of  the  second  order  singularity  present 
in  the  denominator  when  syj=0. 

&{x,  y )  _  f*1 

V  ~  QnJ-i 

(33) 

{/l  "  +/y  7  + 

The  singular  cV|  integration  in  equation  (33)  will  have  to  be  performed 
by  taking  the  finite  part  of  an  infinite  integral  in  the  Mangier  sense 
(32:13)  and  o  in  this  instance  indicates  a  small  distance  away  from  y. 
Note  that  the  Mach  nunber  dependence  has  been  omitted  fran  the  kernel 
function  because  the  flow  conditions  have  been  specified  to  be 
inca.pressibl e .  However,  the  procedure  described  in  this  chapter  is 
applicable  to  compressible  cases  as  well,  provided  the  proper  form  of 
the  kernel  function  is  used.  The  airfoil  surface  is  divided  into  tj- 
integration  regions  to  facilitate  the  integration  process:  the  region 
containing  the  singularity  will  be  referred  to  as  Region  II,  the  region 
to  the  left  of  this  Region  I,  and  to  the  right,  Region  III  as  shown  in 
Figure  6.  As  mentioned  previously,  the  modified  kernel  function 
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Figure  6:  Integration  Regions 

possesses  a  finite  jump  discontinuity  as  x-£-  0  and  y-n-  0.  Otherwise 
the  integration  poses  no  extreme  challenges.  The  wing  is  now  divided 

into  N  panels  as  shown  in  Figure  7.  The  number  of  panels  in  the  x- 
di  recti  on  is  denoted  by  Nj  and  in  the  y -direct  ion  by  Nj  sue*  that  N=Nj  X 
Ny.  The  pressure  coefficient  information  for  the  entire  panel  (i,j) 
will  *je  concentrated  at  the  point  ({ ,tj)  located  at  the  one  quarter  chord 
and  mid-span  position  on  the  panel  and  the  downwash  information 
concentrated  at  the  panel  three  quarter  chord  and  mid-span  point  (x,y) 
as  shown  in  Figure  8.  The  convention  of  utilizing  the  one- quarter/ 


Figure  7:  Discretization  of  Airfoil  Surface 

three-quarter  chord  rule  for  doublet  and  downwash  location  has  been  used 
successfully  for  years  by  both  the  vortex  and  doublet  lattice  methods 
(20:2)  but  the  rigorous  analytical  justification  for  such  an  assumption 
is  lacking  to  date  (20:325-342). 

Hence,  the  pressure  coefficient  distribution  can  be  written  in  the 
following  form  shown  in  equation  (34). 

AP(£,ti;/c)  •  Pij{k)b(t-t1,T\-y)j)  Regional,  III 

(34) 

AP($,T|;k)  =  Pi(ic)«  ($-?i)  Region  II 
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Figure  8:  Airfoil  panel 

Using  equation  (34)  for  the  pressure  coefficient  distribution  and  a 
similar  type  of  fonrulatian  for  the  dawnwash  distribution  produces  the 
following  system  of  linear  equations  relating  the  coefficients  of  the 
downwash  and  pressure  coefficient  distributions. 

In  equation  (35)  the  A^  indicates  the  area  of  the  (i,  j)  panel, 
indicates  the  undetermined  pressure  coefficient  strength  on  the  (i,  j) 
panel  and  is  the  length  of  the  panel  (i,  j)  along  the  spanwise 
coordinate  strip  containing  the  singularity.  The  quantity  w(xD,y1) 
indicates  the  downwash  strength  on  the  (n,m)  panel.  This  approach  still 
requires  the  evaluation  of  the  finite  part  of  an  infinite  integral  but 
it  does  avoid  having  to  treat  the  jurp  discontinuity  in  the  modified 
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kernel  function  because  the  locations  of  the  downwash  and  pressure 
coefficient  sample  points  are  not  coincident  and  hence  Xq*0. 

The  Mangier  integration  can  be  handled  in  one  of  two  ways  as 
mentioned  previously.  The  first  method  by  Ueda  (27:350)  involves 
determining  an  equivalent  non-singular  expression  for  the  kernel 
function  at  y-^=0  based  upon  a  discrete  treatment  of  the  Mangier 
integration  and  an  asynptotic  expansion  for  small  argument  of  the  kernel 
function.  The  second  approach  is  to  utilize  a  modified  Gauss  quadrature 
rule  to  capture  the  finite  part  of  the  integral  (25:1196)  in  which  the 
kernel  function  is  evaluated  at  the  abscissa  points  by  any  accurate 
means.  Both  methods  were  compared  and  produced  comparable  results  using 
as  few  as  six  abscissa  points  with  the  quadrature  rule.  The  former 
method  was  chosen,  however,  because  it  is  more  consistent  with  the 
sampling  of  the  kernel  function  for  the  other  two  regions  of  the 
airfoil.  This  step  collapses  the  present  methodology  into  the  doublet 
point  method  (27:349). 
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The  method  used  in  this  thesis  permits  the  Mangier  integration  to 
be  performed  by  any  appropriate  means.  Another  distinction  from  the 
doublet  point  method  (27:349)  is  the  manner  in  which  the  kernel  function 
is  evaluated.  The  present  method  computed  the  kernel  function  in  the 
non-singular  regions  by  first  computing  the  modified  kernel  function, 
equation  (5)  by  methods  described  in  Section  II  and  then  dividing  by  the 
double-pole  singularity  term  (equation  (4)).  The  kernel  function  in  the 
singular  region  was  handled  by  the  approximation  given  in  equation  (46). 
The  doublet  point  method  (27:349)  utilizes  the  asymptotic  expansions, 
equations  (39)  through  (44)  to  compute  the  kernel  function  throughout 
the  domain  of  the  airfoil .  The  series  converge  at  different  rates 
depending  on  the  values  of  the  arguments,  requiring  special  treatment  of 
summation  computationally.  Rather  than  work  with  this,  the  alternate 
method  was  used  which  did  not  require  any  series  evaluations. 

Ueda  has  shown  (26:169)  the  kernel  function  can  be  written  as 
followi  for  incompressible  flow  conditions 

K[k,x-t,  s(y-ti)  ]  -  e'ixx°B(k,x-^,  s{y-i\) )  (36* 
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where  the  function  8(k,x-£ ,s(y-^ ) )  is  defined 


B(k,x-Z,  s(y-r\) ) 


r*-t _ e1** _ 

J—  [va+sa  (y-t)) 3  ]3/a 


(37) 


If  this  integral  is  separated  into  its  real  and  imaginary  parts,  the 
real  part  contains  a  double-pole  singularity  in  ^  as  n-  0  as  shown  in 
equation  (38)  which  was  obtained  by  making  a  change  of  variables  on 
equation  (37). 


Bik.XQ,  8  (y- n) ) 
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Ueda  has  expanded  both  the  real  and  imaginary  parts  of  equation  in  an 
asymptotic  series  about  small  ks(y-^).  The  series  is  written  as  the 
following  (21:170) 
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One  attractive  feature  of  these  series  is  that  every  singular  part  is 
explicitly  written  as  an  initial  term  (26:169).  These  expressions  can 
be  inserted  in  the  q  integral  and  the  integrations  can  be  performed  in 
the  Mangier  sense. 


f'BK{k,x0lsy0)  d(sy0)  -  ^[/"^  +  j* B^k,^,  By0)  d{By0)  -  ] 

*o>0 

(45) 

Ueda  has  shown  an  equivalent  expression  for  Bg  after  including  the 
effects  of  the  Mangier  integral  performed  in  a  discrete  sense  can  be 
written  as  the  following  (27:350) 


BK(k,XQ,sy0)  -  Bn{k,  -*0,0)  -  +  y  -  -|| 


(46) 


(sy0<o,  Xg>0) 
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This  provides  a  method  of  evaluating  the  kernel  function  when  syg=0  and 
can  be  substituted  into  the  second  summation  expression  in  equation 
(35).  Replace  the  integration  in  the  second  expression  in  equation  (35) 
with  a  discrete  sampling  of  the  kernel  function  as  done  previously  in 
equation  (34).  This  permits  the  panel  numbering  system  (i,j)  to 
collapse  to  a  single  because  it  is  no  longer  required  to  track  specially 
the  spanwise  coordinate  of  the  modified  kernel  function  for  integration 
purposes  as  there  exists  a  way  to  evaluate  the  kernel  function  at  sy^=0 
and  no  special  integration  treatment  is  now  required  for  the  kernel 
function.  Hence,  the  (i, j)  panel  will  be  relabeled  as  the  i-th  panel 
with  a  doublet  point  at  and  downwash  at  point  x^,y^.  This 
sinplification  allows  the  system  of  equations  in  equation  (35)  to  be 
expressed  in  a  much  simpler  form 

*  K&  (47> 

&  -{Py} 

*  -  {**}  ■  (49) 

(50) 

where  i  and  j  can  index  from  1  to  N=Nj  X  Ny  panels.  The  ->  indicates  a 
vector  quantity  and  K  distinguishes  the  matrix  quantity  from  the 
continuous  quantity.  This  system  of  equations  (47)  through  (50)  is 
presented  as  the  doublet  point  method  (27:349).  This  linear  system  of 


51 


equations  can  be  solved  given  an  assured  downwash  modeshape  vector  . 
The  beauty  in  this  method  is  evidenced  by  not  having  to  presune  an 
appropriate  form  for  the  pressure  coefficient  distribution.  This  method 
relies  on  the  mathematics  embodied  inside  the  kernel  function  to 
determine  the  correct  pressure  coefficient  distribution.  An 
illustration  of  the  ability  of  the  procedure  to  compute  the  pressure 
coefficient  distribution  is  seen  in  Figure  9  below.  Most  methods 
discussed  previously  such  as  the  vortex  and  doublet  lattice  methods  and 
all  assuned  mode  methods,  involve  assuring  some  form  for  the  pressure 
coefficient  distribution.  This  distribution  is  generally  singular, 
creating  even  more  difficulties  using  one  of  these  other  approximate 
methods.  The  pressure  coefficient  Ap  shown  is  for  a  thin  rectangular  of 
aspect  ratio  two  in  steady  flow.  Hie  pressure  distribution  on  an 
airfoil  of  this  type  under  steady  flow  conditions  should  possess  two 
distinct  features:  first,  the  flow  should  exhibit  a  singularity  of  the 
type  e'l,  6-.  0  at  the  leading  edge  of  the  airfoil  (32:5)  and  the  flow  to 
go  to  zero  at  the  wing  tips  with  an  infinite  slope  as  el  6-  0  (32:5). 
The  first  type  of  behavior  is  demonstrated  in  Figure  9,  becoming  more 
pronounced  as  the  number  of  chordwise  panels  Nj  is  increased.  The 
second  type  of  behavior  is  seen  in  Figure  10  where  Cj  4  is  shown  as  a 
function  of  the  nuiber  of  spanwise  panels,  Nj.  The  coefficient  of  lift 
is  calculated  using  equation  (60)  which  will  be  introduced  shortly.  The 
behavior  is  more  pronounced  as  the  ntvrber  of  panels  is  increased,  as 
observed  in  the  previous  figure. 
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CHORDWISE  POSITION 


Figure  9:  Dependence  on  Nj 

The  total  pitching  lift  coefficient  for  the  case  Ny=10  is  2.78 
compared  with  the  previously  reported  value  of  2.65  when  Ny=20.  This 
demonstrates  the  dependence  between  the  total  pitching  lift  coefficient 
and  the  number  of  panels  in  the  spanwise  direction.  The  total  pitching 
lift  coefficient  for  the  Nx=5  case  was  2.65  ccnpared  to  the  1^=10  case 
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with  a  value  of  2.64.  This  shows  the  total  pitching  lift  coefficient  is 
little  changed  by  the  number  of  panels. 


SPANWISE  POSITION 


Figure  10:  Dependence  on  Nj 

These  exanples  serve  to  show  the  increased  accuracy  in  the  local 
pitching  lift  coefficient  achieved  by  the  larger  number  of  panels  in  the 
chordwise  direction  and  little  effect  the  number  of  chordwise  panels  has 
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cxi  the  same.  The  unsteady  cases  also  have  a  similar  type  of  dependence 
upon  the  number  of  s panwise  and  chordwise  panels. 

This  procedure  is  not  compared  to  actual  test  data  here,  but  is 
compared  to  other  approximate  methods.  Test  data  such  as  that  in 
reference  (34)  is  difficult  to  carpare  with  at  lower  reduced 
frequencies.  This  is  due  to  the  very  high  velocities  necessary  to 
achieve  low  reduced  frequencies  (k=#t^/V)  inside  the  wind  tunnel.  These 
higher  velocities  result  in  high  mach  numbers,  violating  the 
incompressible  flow  assumptions  upon  which  the  procedure  is  based.  Even 
Mach  number  corrections  are  not  very  helpful  because  the  higher  Mach 
numbers  feature  shock  waves  traveling  along  the  airfoil  which  is  a  non¬ 
linear  problem  not  handled  by  this  theory.  Comparisons  of  the  computed 
steady  values  of  total  pitching  lift  coefficient  with  those  of  modified 
lifting  line  theory  (20:338)  and  an  approximate  method  by  Graham  (7:93) 
show  a  maximum  error  of  approximately  six  percent  and  a  minimum  error  of 
about  0.3  percent  for  aspect  ratios  ranging  from  two  through  ten  with 
Nj=5  and  Nj=20.  Although  these  appear  large,  the  resul ts  are 
respectable  for  unsteady  aerodynamic  analysis.  The  unsteady  results 
over  a  range  of  reduced  frequency  from  zero  to  one  for  the  aspect  ratio 
two  cases  were  compared  with  those  of  another  approximate  method  by 
Lawrence  (15:769-781)  and  the  agreement  was  good. 

There  is  a  desire  to  use  a  large  number  of  elements  to  achieve 
higher  accuracy.  However,  this  increases  the  size  of  the  linear  system 
and  significantly  increases  the  solution  time,  especially  when  the 
computations  must  be  done  over  a  wide  range  of  reduced  frequencies.  The 
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aspect  ratio  two  case  was  solved  using  a  100  panel  (5  X  20)  wing  and 
aspect  ratio  ten  was  solved  using  a  600  panel  (5  X  100)  airfoil.  The 
dependence  upon  the  Nj  for  each  aspect  ratio  was  determined  by  keeping 
the  Nj/AR  ratio  constant.  This  method  presented  in  this  section  and  the 
doublet  point  method  (27:349)  is  much  less  time  consuming 
computationally  than  any  of  the  procedures  mentioned  earlier.  Another 
option  that  exists  for  wings  with  s panwise  symmetry  is  to  reduce  the 
number  of  degrees  of  freedom  by  utilizing  a  constraint  matrix  (27:351) 

If  the  pressure  coefficient  distribution  can  be  represented  by  a  smaller 
number  of  degrees  of  freedom 

&  -  Qp 

where  Q  is  the  constraint  matrix.  This  allows  equation  (47)  to  be 
written  as 

Qt#  *  [QTK(j\p  <52> 

where  the  number  of  degrees  of  freedom  is  reduced  by  an  amount 
determined  by  the  Q  matrix.  Generally,  for  rectangular,  spanwise 
symmetry,  the  number  of  degrees  of  freedom  will  be  reduced  by  a  factor 
of  one-half. 

An  important  feature  of  this  procedure  is  direct  input  into  finite 
element  types  of  aeroelastic  analyses.  Once  the  structural  grid  is 
established,  this  routine  can  be  used  to  obtain  reasonably  accurate 
structural  applied  loads  by  computing  the  pressure  coefficient  vector 
for  the  structure  and  multiplying  by  ipV2S.  The  finite  element 
technique  can  then  be  used  to  solve  the  aeroelastic  problem. 
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Ashley  and  Boyd  (2:8)  have  demonstrated  the  expression  for  the 
modified  kernel  function  in  equation  (5)  is  an  appropriate  form  for 
computing  values  throughout  the  complex  s -plane  provided  special 
precautions  are  taken  in  computing  the  modified  kernel  function  in 
certain  regions  of  the  s-plane.  Hence,  the  generalized  forces  computed 
using  this  method  can  be  treated  as  Laplace  transformed  aerodynamic 
forces.  Modelled  with  fractional  calculus,  the  generalized  forces  can 
now  be  used  more  easily  to  determine  the  forced  response  of  the  airfoil 
which  will  be  discussed  in  Section  VI. 
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VI.  Equivalent  Theodor  sen  Function  for  Three-dimensional  Aerodynamics 

Theodorsen's  function  is  the  transfer  function  relating  the 
unsteady  circulatory  lift  to  the  downwash  of  an  infinitely  thin, 
harmonically  oscillating  plate  in  inviscid,  incompressible,  two- 
dimensional  flow  (25:418).  The  Theodorsen  function  is  a  ratio  of  Bessel 
functions  with  no  known  analytic  inverse  transform  whose  arguments  are 
defined  in  terms  of  the  reduced  frequency  and  taken  as  purely  imaginary 
for  stability  analysis  (4:5).  Swinney  (24:5)  has  developed  a  fractional 
calculus  model  of  the  Theodorsen  function  which  is  valid  in  the  entire 
s -plane  with  the  appropriate  substitution  of  the  dimensionless  Laplace 
variable  s  for  ik.  Unfortunately,  no  direct  functional  analog  to  the 
Theodorsen  function  in  three-dimensional  aerodynamics  has  been 
established.  The  integral  equation  (3)  shows  the  transfer  function  is 
part  of  a  spatial  convolution  with  the  transform  of  the  pressure 
coefficient  distribution  on  the  airfoil.  This  transfer  function  is  not 
as  mathematically  clean  as  the  Theodorsen  function  and  little 
investigation  and  progress  has  been  made  toward  the  simplification  of 
the  three-dimensional  unsteady  aerodynamic  transfer  function. 

This  chapter  will  describe  an  approach  taken  to  simplify  the 
determination  of  the  transfer  function  in  three-dimensional 
aerodynamics.  As  stated  above,  this  transfer  function  is  embedded  in 
the  spatial  convolution  in  the  integral  equation  (3).  The  approach 
taken  is  to  express  the  generalized  aerodynamic  forces  in  three- 
dimensional  flow  in  terms  of  the  forces  found  in  two-dimensional  flow  by 
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appropriate  use  of  a  new  function,  the  equivalent  Theodorsen  function 
for  three-dimensional  unsteady  aerodynamics.  The  equivalent  Theodorsen 
function  will  be  developed  and  its  properties  explained.  Additionally, 
since  the  equivalent  Theodorsen  function  exists  in  the  Laplace  domain, 
the  acccnpanying  time  domain  functions  used  in  aeroelastic  analysis  will 
be  developed  as  well. 

The  idea  of  relating  the  two-  and  three-dimensional  aerodynamic 
forces  is  not  new.  Many  different  theories  have  been  authored 
describing  how  to  establish  such  a  relationship.  Watkins,  Runyan  and 
Wool stan  established  that  the  kernel  function  does  collapse  to  the  two- 
dimensional  case  as  the  span  is  increased  to  infinity  (31:713),  which 
may  suggest  the  existence  of  same  type  of  mathematical  link  between  the 
finite  span  wing  and  infinite  airfoil  unsteady  aerodynamic  phenomena. 
Reissner  (21:32)  has  established  an  integral  correction  factor  based 
upon  an  strip  formula  type  of  approach  similar  to  that  of  lifting-line 
theory  which  modifies  the  two-dimensional  Theodorsen  function  to  produce 
the  generalized  forces  on  a  finite  wing.  Another  method  has  been 
developed  by  Lawrence  (15:771)  which  modifies  the  Theodorsen  function  to 
produce  three-dimensional  effects.  These  methods  are  based  upon 
different  underlying  assumptions  regarding  the  connection  between  the 
two- and  three-dimensional  flow  conditions.  Reissner 's  theory 
establishes  the  non-circulatory  portion  of  the  lift  on  aui  infinite  wing 
remains  unchanged  as  the  span  becomes  finite  (6:389).  The  non- 
circulatory  flow  produces  no  wake  and  therefore,  the  induction  effects 
associated  with  it  are  relatively  small.  In  regions  of  the  wing  where 
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the  induction  effects  are  small  compared  to  the  purely  circulatory  ones, 
this  is  certainly  true.  However,  this  is  not  true  in  the  immediate  wing 
tip  region  because  of  the  larger  induction  effects  from  the  tip  vortices 
present.  Reissner  has  also  established  the  uninportance  of  induction 
effects  on  the  accurate  prediction  of  airloads  on  wings  for  large  values 
of  reduced  frequency  since  the  nan-circulatory  effects  are  proportional 
to  k2  (6:389).  However,  induction  effects  are  dominant  for  small  values 
of  k  and  for  small  aspect  ratio  wings.  These  observations  and  theories 
led  to  the  postulation  of  the  equivalent  three-dimensional  Theodorsen 
function. 

Development 

Two-dimensional  unsteady  aerodynamic  equations  are  an  appropriate 
starting  point  to  begin  a  search  for  an  appropriate  form  of  the 
equivalent  Theodorsen  function.  Lift  generated  by  pure  pitching  motion 
was  selected  as  a  generalized  aerodynamic  force  with  which  to  compare 
the  two-  and  three-dimensional  unsteady  aerodynamics.  The  lift  per  unit 
span  of  an  airfoil  in  pure  pitch  is  given  by  (6:272) 

L  -  npi?o  [Vd]  +  2xpUb0C  (k)  [V«  ♦  -^4]  (53) 

where  C(k)  is  the  Theodorsen  function  for  two-dimensional  aerodynamics. 
The  first  term  in  equation  (53)  represents  the  non- circulatory  portion 
and  the  second  term  represent  the  circulatory  portion  of  lift  per  unit 
span.  In  this  case  the  dcwnwash  function  w(x,y,t)=a(x,y,t)  has  taken  on 
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the  following  form  to  insure  positive  angle  at  the  leading  edge  of  the 
airfoil  and  to  be  consistent  with  equation  (53) 

w(x,y,t )  -#  (x,y,  t)  -  -a0xelvc  (54) 

Noting  that  lift  per  unit  span  for  the  finite  rectangular  airfoil  can  be 
approximated  by  the  following  expression  while  still  being  consistent 
with  the  original  assumptions  of  the  integral  equation  (1),  it  is 
possible  to  directly  relate  the  integral  equation  for  three-dimensional 
flow  to  the  Theodorsen  function. 


L  -  -|pV*C^«(2i>0) 

a  _  *(x,y,  t) 

V 


(55) 


Here  the  T(k)  has  been  substituted  for  the  usual  C(k)  and  will  now  be 
referred  to  as  the  equivalent  Theodorsen  function  for  three-dimensional 
unsteady  aerodynamics  because  the  coefficient  of  lift,  C^,  is  computed 
using  the  three-dimensional  theory.  Substituting  these  expressions  into 
equation  (53)  and  simplifying  yields 


-  *  ji/e  +  arwji  +  (56) 

No. e  that  this  model  presunes  the  changes  between  two-and  three- 
dirensional  flow  exists  purely  in  the  equivalent  Theodorsen  function. 

Thf.  initial  model  was  constructed  with  no  constraints  placed  on  the  ncn- 
circulatory  and  circulatory  lift,  alleging  the  two  to  vary  independently 
of  one  another.  The  initial  model  was  found  to  be  suboptimal .  Next, 
the  model  in  equation  (56)  in  which  the  ncn-circul atony  lift  is  fixed 
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and  the  circulatory  lift  is  allowed  to  vary  was  used.  This  model  was 
shown  to  yield  satisfactory  results  and  is  an  incomplete  confirmation  of 
Reissner’s  conjecture  mentioned  earlier  which  stated  that  the  nan- 
circulatory  lift  does  not  change.  This  will  elaborated  upon  further 
after  the  fractional  calculus  model  of  the  equivalent  Theodorsen 
function  is  presented  in  the  next  section.  Given  the  relationship  shown 
in  equation  (56),  discrete  values  of  T(k)  can  be  ccnputed  from  the 
discrete  values  of  the  coefficient  of  lift  and  modelled  using  fractional 
calculus.  Given  a  convenient  form  for  T(k),  time  domain  response  will 
be  calculable  as  is  the  case  with  C(k).  It  should  be  noted  that  pure 
pitching  motion  is  only  one  form  of  motion  for  which  this  relationship 
should  exist.  Another  form  being  pure  plunge,  for  example. 

Fractional  Calculus  Model  and  Results 

A  fractional  calculus  model  for  the  two-dimensional  Theodorsen 
function  (24:5)  suggested  a  fractional  calculus  form  for  the  equivalent 
three-dimensional  function  since  the  two  have  been  presumed  to  be 
mathematically  related  as  shown  in  equation  (56).  Also,  as  shown 
earlier,  the  kernel  function  for  incompressible  flow  conditions  was 
shown  to  possess  fractional  derivative  properties  which  further 
suggested  the  resulting  lift  and  hence,  the  equivalent  Theodorsen 
function  T(k),  should  have  these  same  trends. 

The  equivalent  Theodorsen  function  is  highly  nonlinear  and  there 
is  a  number  of  possible  forms  of  models.  In  this  case,  Swinney’s  model 
of  the  Theodorsen  function  was  used  as  a  suitable  form  with  which  to 
begin.  Before  the  modelling  can  begin,  the  coefficient  of  lift  must  be 
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computed.  The  next  paragraph  will  elaborate  on  the  exact  method  used  to 
determine  the  coefficient  of  lift. 

The  mode  shape  for  pure  pitching  about  midchord  can  be  written  in 
a  form  consistent  with  equation  (19)  as 

a(x.y.t)  -  <57> 

bo 

where  OgX  is  the  shape  function  describing  the  position  of  all  points  on 
the  airfoil  experiencing  pure  pitching  motion  and  the  minus  sign  is 
included  to  produce  positive  angle  for  pos'^ive  Og  at  the  leading  edge. 
The  subscript  a  makes  reference  to  the  pitching  amplitude  function 
q j  ( t ) .  There  is  no  y  dependence  in  the  shape  function  because  there  is 
pure  rigid  body  pitching  and  no  spanwise  bending  or  torsion  is  present. 
This  downwash  can  be  cast  into  the  form  in  equation  (49)  as 

*  -  0  -  -  {-«o(1+i**i)}  (58) 

Using  the  relationship  defined  previously  in  equation  (50)  for  K,  the 
unknown  pressure  coefficient  can  now  be  solved  via  standard  linear 
system  procedures.  The  K  matrix  is  well  conditioned  and  solutions  were 
obtained  using  the  IMSL  routine  LSAOG  which  solves  complex  general 
systems  of  linear  equations  by  Gaussian  elimination  with  iterative 
refinement  (9:31).  Given  the  pressure  coefficient,  the  following 
definitions  can  be  used  to  obtain  the  ccnplex  lift  coefficient  per  unit 
span  and  total  lift  coefficient  as  shown  in  equations  (59)  and  (60).  The 
coefficients  of  lift  in  equations  (59)  and  (60)  above  can  be  treated  as 


(59) 


ci  (yj) 


PiJ&i 


c  =  1  y  ( 


ci  (y^)  ^ 


(60) 


derivatives  with  respect  to  angle  of  attack,  a,  because  the  coefficients 
of  pressure  upon  which  they  are  based  are  only  known  to  an 
arbitrary  magnitude,  Oj,  of  the  pitching  motion. 

The  equivalent  Theodor sen  function  modelling  was  acccnplished 
using  the  same  procedure  as  was  performed  for  modelling  the  modified 
kernel  function,  the  only  difference  being  the  real  part  was  used  for 
regression  in  this  instance  versus  the  imaginary  part  in  the  prior  case. 
The  real  part  inexplicably  seemed  to  produce  more  consistent  results  in 
this  case.  The  functional  form  used  to  model  the  equivalent  Theodor sen 
function  is 


T(k.AR) 


T  +  *  (AS)  (i*)M*° 

0  1  +  b(AR)  (i*)'*<") 


-  gr(AR)  (Ik)  +  f(AR)  (Ik) 2 


(61) 


where  Tq  is  given  by  the  steady  value  of  #  as  determined  in  equation 
(56)  divided  by  %.  The  results  once  again  yielded  a  family  of  curves, 
the  parameters  in  this  instance  being  functions  of  aspect  ratio  only. 
The  accuracy  of  the  modelling  of  the  equivalent  Theodorsen  function  is 
shown  below  and  in  Appendix  D.  The  equivalent  Theodorsen  function 
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REDUCED  FREQUENCY 


Figure  12:  Imaginary  Parts  of  Theodorsen  Functions 

Theodorsen  function  T(k)  tends  toward  the  two-dimensional  C(k)  in  the 
limit  as  the  aspect  ratio  is  increased  from  two  to  ten;  this  was  not 
verified  or  proven  in  this  work. 

T(ik)  =  r0  +  — .  _  gr(ijt)  +  f(ik)3  (6 

0  1  +  b(lk)* 


U=0.904 
f =0.085 
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b=3.32 
g=0.017 


-  Model  of  Ct,oipho 
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REDUCED  FREQUENCY 


Figure  13:  for  Wing  of  Aspect  Ratio  10 


It  would  be  interesting  to  see  the  results  for  higher  aspect  ratio  such 
as  twenty  or  fifty.  This  type  of  analysis  was  not  pursued  further  here 
because  of  a  lade  of  available  ccnputer  memory.  This  limiting 
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relationship  was  not  proven  here  but  does  suggest  a  topic  for  further 
study.  As  mentioned  earlier,  Reissner  (22:14-18)  has  developed  a 
relationship  between  the  two-dimensional  Theodors en  function  and  three- 
dimensional  unsteady  aerodynamics.  This  relationship  consists  of 
modifying  the  Theodor sen  function  by  multiplication  with  a  factor 
consisting  of  an  integral  factor  and  a  combination  of  Bessel  functions. 
The  integrals  in  Reissner' s  relationship  are  quite  difficult  to  evaluate 
and  have  been  tabulated  for  easier  use  (22:70-72).  It  is  believed  the 
behavior  within  the  fractional  calculus  model  may  fully  capture  the 
adjustment  Reissner  made  to  the  Theodorsen  function  with  his  correction 
factor.  The  theory  by  Reissner  presumes  the  non-ci rcul atory  effects  of 
the  flow  remain  unchanged  between  the  two-  and  three-dimensional  flow 
(6:389).  The  mathematics  of  the  modelling  process  suggested  this 
phenomenon  as  well.  In  numerous  attempts  to  force  the  three-dimensional 
effects  to  be  contained  in  an  additional  non-ci  rcul  atory  term,  the 
models  would  not  converge  uniformly  for  both  the  real  and  imaginary 
parts  of  the  coefficients  of  lift.  Many  different  types  of  additive 
non-ci rcul atory  terms  and  many  different  trial  vectors  were  tried  with 
no  success .  It  was  only  when  the  flow  modifications  were  restricted  to 
the  circulatory  portion  of  the  flow  that  good  convergence  for  both  the 
real  and  imaginary  parts  was  obtained.  The  Theodorsen  function  is  the 
transfer  function  for  the  circulatory  portion  of  two-dimensional  flow. 
The  relative  contribution  of  circulatory  lift  to  the  coefficient  of 
lift,  the  quantity  [l+ik/2]  in  equation  (53),  which  multiplies  the 
Theodorsen  function  for  pure  pitching  motion  as  shown  in  equation 
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remained  unchanged  and  all  differences  between  the  two-  and  three- 
dimensional  flow  were  placed  in  the  equivalent  Theodorsen  function. 

Another  interesting  result  from  this  analysis  is  the  relationship 
established  among  the  parameters  of  the  model.  This  relationship  can  be 
used  to  determine  other  values  of  the  equivalent  Theodorsen  function  for 
rectangular  wings  of  aspect  ratios  between  two  and  ten.  The  curves 
relating  the  parameters  are  contained  in  Appendix  E  with  a  sample  shown 
in  Figure  14  and  Figure  15.  These  results  are  good  only  for  rectangular 
planforms  but  it  is  believed  such  a  relationship  might  exist  tor  other 
planforms  such  as  tapered  and  swept -tapered. 

The  parameter  a  was  well  behaved  but  the  y  parameter  was  more 
complicated.  More  is  to  be  gained  mathematically  by  viewing  the 
polynomial  model  of  these  parameters  as  the  first  portion  of  an  infinite 
power  series  of  some  unknown  function  which  fully  describes  the 
behavior.  A  cubic  spline  is  plotted  through  the  y  points  as  the  order 
of  the  polynomial  was  continually  increasing.  The  coefficients  of  the 
polynomial  were  decreasing  toward  zero  quite  rapidly,  typical  of  a 
convergent  series.  It  is  postulated  that  the  unusual  behavior  of  the 
parameter  y  may  smooth  out  if  the  modelled  range  of  reduced  frequency  is 
increased. 
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Figure  14:  Parameter  a  vs.  Aspect  Ratio 


MODEL  PARAMETER  Mu(AR) 


Figure  15:  Parameter  u  vs.  Aspect  Ratio 


VII  •  Time  Dcrrain  Responses  Based  Upon  Equivalent  Theodorsen  Function 


In  general,  it  is  difficult  to  get  an  accurate  time  domain 
representation  of  unsteady  aerodynamic  loads  using  the  three-dimensional 
pressure-downwash  integral  equation  (1).  Usually,  the  solution  of  this 
equation  is  only  presented  for  frequency  only  since  most  of  the  work 
with  the  kernel  function  has  been  restricted  to  stability  or  flutter 
analyses  which  are  only  defined  along  the  positive  imaginary  axis.  In 
two-dimensional  unsteady  aerodynamics,  the  Theodorsen  function  has  been 
approximated  with  a  model  that  has  an  elementary  inverse  Laplace 
transform  to  permit  forced  response  analysis.  As  this  is  the  first  time 
the  circulatory  properties  of  the  three-dimensional  aerodynamics  have 
been  cast  into  a  mathematically  tractable  form,  this  will  be  the  first 
time  the  forced  responses  for  wings  of  finite  span  will  be  written  in  a 
mathematically  accurate  and  compact  form. 

Wagner  Function 

The  advantage  of  this  fractional  calculus  model  is  the  reasonably 
compact  time  domain  representation  of  the  unsteady  aerodynamic  loads. 

The  other  approximate  methods  used  to  compute  the  generalized 
aerodynamic  forces,  such  as  Reissner's  (21:1-39)  and  Lawrence's  (15:769- 
781),  do  not  have  this  unique  and  beneficial  feature.  The  time- 
dependent  lift  is  posed  in  terms  of  the  Wagner  function  (6:285),  $(t), 
which  is  defined  in  equation  (63).  The  Wagner  function  represents  the 
lift  resulting  from  a  unit  step  change  in  angle  of  attack  (sometimes 
called  indicial  lift). 
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(63) 


♦  <t) 


= 


'  C(s)  • 
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The  fractional  calculus  model  of  the  equivalent  Theodorsen  function, 
T(k),  in  equation  (61)  can  be  represented  in  the  Laplace  domain  as 


T(s)  *  T0  +  — -  —  -  g(a)  +  f(s) 2  (64) 

0  1  +  b(s)* 


where  the  dimensionless  Laplace  variable  s=§fc^/V=ik  indicates  the  change 
into  the  Laplace  plane,  s  the  general  Laplace  variable,  and  where  the  AR 
dependence  of  the  parameters  is  not  shown  for  compactness.  Equation 
(64)  can  manipulated  into  a  form  more  suitable  for  inverse  Laplace 
transformation 


T{b) 


t0 


a  +  a  1 
b  b  b(s)*  +  l 


g(a)  +  f{s)2 


(65) 


By  recognizing  the  binomial  series  present,  this  can  be  further 
expressed  by  the  following  power  series 

r<5)  tf<s,a  (6 

These  simplifications  lead  to  the  following  Laplace  domain  expression 
for  the  Wagner  function 


r(a) 

$ 


To 
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_  a  A  /-i\°  1  1 
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(67) 


Making  a  substitution  for  s  the  Wagner  function  can  be  written  in  terms 
of  the  Laplace  variable  S 
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Taking  the  inverse  Laplace  transform  yields 


♦  (t)  =  T0u.x(t) 


A  T  l -ir(  vtY°  1 

b  /rj  l  b  I  {  b0  j  r(nj2  +  i) 
-  fu+1(t) 


(69) 


The  term  u^  is  a  unit  step  function,  the  Ug  term  is  the  unit  delta 
function,  and  the  u^  term  is  the  unit  doublet  function.  The  doublet 
function  is  the  derivative  of  the  delta  function.  The  p  order  Mittag- 
Leffler  function  is  defined  as  (15:102) 


x° 

r(l+pn) 


(70) 


The  p  order  Mittag-Leff ler  function  shown  in  equation  (70)  can  be 
considered  as  the  generalized  p  order  exponential  function  (4:12)  as 
setting  p  equal  to  one  produces  the  Taylor  series  for  the  exponential 
function  e1.  Using  this  notation  to  represent  the  time  series  in  the 
Wagner  function  produces 


♦  (fc)  =  T0u_x(t)  - 


bE» 


-i 

b 


(«rl  -  * 

•(*)' 


9u0  ( t) 


(71) 


fu+1(t) 


74 


The  Wagner  function  is  shown  in  terms  of  dimensionless  time,  t=Vt/t^,  in 
equation  (72). 


-  ^sruo(t) 

(72) 

♦  (if)  fu.i<F> 

This  expression  is  similar  to  the  Wagner  function  determined  by  Swinney 
(24:22)  with  the  addition  of  the  discontinuous  Dirac  delta  function  and 
doublet  function.  The  approximation  to  Wagner's  function  in  three- 
dimensional  flow  (72)  is  shown  with  the  Wagner  function  model  developed 
by  Swinney  (24:22)  for  two-dimensional  flow  in  Figure  16.  It  should  be 
noted  that  only  the  continuous  portions  of  the  Wagner  function 
representation  in  equation  (72)  are  shown  in  Figure  16.  The  figure 
illustrates  the  response  to  a  step  input  in  angle  of  attack  reaches  a 
steady-state  value  faster  in  three-dimensional  flow  than  in  two- 
dimensional  flow.  This  is  expected  as  the  finite  wing  has  wing  tip 
effects  working  to  dampen  motion  in  addition  to  the  trailing  edge 
effects  it  has  in  curnnon  with  the  infinite  wing. 


4>(t)  -  T0u_x(t)  -  | 


-  C 
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WAGNER  FUNCTION 


Time  Dependent  Lift  in  Terns  of  Wagner  Function 

The  lift  for  a  wing  can  be  developed  using  the  Wagner  function 
quite  easily.  The  lift  per  unit  span  for  an  airfoil  was  given  in 
equation  (53)  for  pure  pitching  motion.  Generalizing  this  equation  to 
include  vertical  translation  as  well  and  constructing  the  Laplace 
transform  yields  equation  (73). 

L(§)  =  *p£>oa$£[J5(t)  +  VA(t)] 

(?3) 

+  2*pVi>0S£  A(t)  ♦  Va  +  -^d(t)  T(s) 

mk 

Substituting  the  expression  in  equation  (64)  for  T(s)  and  taking  the 
inverse  Laplace  transform  produces  equation  (74). 

L(t)  -  xpJ?|[^(t)  +  Vd(t)] 

+  2%pVb0  |t0  -  -|J  h(C)  +  Vt(t)  +  ^yd(t) 

-  2xpJ fi(c-v)  +  Vd(t-x)  +  ^~t(t-x)  dx  (74) 

Jo  ax  2 

♦  2y°  fft-^rA(t-x)  ♦  v*{t-x)  *  ( t-t)  ldt 

V  Jo  dx*[  2  j 

-  2*pVb0fJ  h(x)  *  VA(x)  *  y*(T)  ^[^(  —  ^"T))l>  dx 

The  dot  over  the  Mittag-Leffler  function  implies  the  continuous  portion 
of  the  derivative  as  the  discontinuous  portions  are  contained  in  the 
second  term. 
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Kiissner  Function 

A  similar  time-dependent  lift  function  is  the  Kiissner  function 
which  describes  the  lift  created  by  a  sharp-edged  gust  striking  the 
leading  edge  of  an  airfoil.  The  Kiissner  function  $(t)  is  defined  as 
(6:287) 

♦  <t)  -  S|-1{e-*[(J0<a)  -I1{S))T{a)  +2^5)]}  (75) 

Taking  the  inverse  transform  of  equation  (75)  yields  the  time  domain 
representation  of  the  Kiissner  function. 

Hence,  for  rectangular  geometry  wings  between  aspect  ratio  two  and 
ten,  it  is  beneficial  to  utilize  the  equivalent  Theodorsen  function 
inserted  into  equations  of  motion  developed  for  two-dimensional  flow  to 
compute  aerodynamic  loads  and  responses  for  three-dimensional  flow 
conditions.  With  the  development  of  the  equivalent  Theodorsen  function 
for  finite  span  wings,  the  three-dimensional  unsteady  aerodynamic 
analysis  of  response  to  arbitrary  motion  can  be  easily  obtained,  a  task 
which  was  extremely  difficult  to  accomplish  with  the  method.-  described 
previously.  The  advantage  of  the  equivalent  Theodorsen  function  is  that 
it  permits  the  aerodynamic  loads  to  be  written  in  a  mathematically 
compact  form  and  convenient,  closed  form  inverse  Laplace  transforms 
exist  for  the  equivalent  Theodorsen  function,  making  its  use  quite 
attractive. 
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VIII .  Conclusions 


This  thesis  has  investigated  two  alternative  approaches  to 
predicting  the  aerodynamic  loads  produced  by  time-dependent  motions  of 
thin  wings  in  rectilinear  subsonic  flight.  The  first  method  was  an 
attempt  to  directly  integrate  the  pressure-down wash  integral  equation  in 
three-dimensional  unsteady  aerodynamics.  This  method  consisted  of 
finding  a  model  of  the  three-dimensional  kernel  function  to  capture  the 
behavior  in  a  more  mathematically  convenient  form  to  facilitate  the 
spatial  convolution.  The  model  was  sought  in  a  fractional  calculus  form 
to  capture  better  the  frequency  dependent  properties  and  to  make  the 
analytic  continuation  from  the  imaginary  axis  into  the  entire  complex 
s -plane  easier.  It  was  found  that  the  canpressibility  effects  within 
the  kernel  function  made  the  modelling  at  higher  Mach  ntmbers  extremely 
difficult.  The  kernel  function  in  incompressible  flow  was  shown  to  be 
modelled  by  a  simple  four-parameter  model  throughout  most  of  the  domain 
of  small  aspect  ratio  airfoils;  however,  the  small  and  large  asymptotes 
of  reduced  frequency  have  not  yet  been  verified.  Initial  investigations 
suggest  this  should  be  successful  as  well.  Additionally,  the 
determination  of  appropriate  behavior  for  the  model  in  the  complex  s- 
plane  has  not  been  verified  yet  either  should  be  accomplished  in  the 
future.  The  four-parameter  model  began  to  worsen  as  the  values  of  the 
s panwise  variable  were  increased  above  those  of  aspect  ratio  one-half. 
The  higher  cases  of  aspect  ratio  need  further  investigation  with 
different  forms  of  fractional  calculus  models.  The  inplementation  of  a 
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complex,  nonlinear,  least-squares  regression  algorithm  with  adaptive 
search  methods  would  greatly  enhance  the  capability  of  examining  these 
models.  The  three-dimensional  kernel  function  contains  fractional 
calculus  behavior  and  it  is  only  a  matter  of  time  before  a  suitable 
model,  mathematically  simpler  in  form  than  the  kernel  function  itself, 
will  be  found  and  successfully  applied  to  the  active  control  system 
design  of  aircraft. 

The  second  approach  taken  to  the  prediction  of  unsteady 
aerodynamic  loads  on  finite  span  rectangular  airfoils  was  to  develop  an 
equivalent  Theodorsen  function.  The  equivalent  Theodorsen  can  then  be 
substituted  into  the  two-dimensional  air  load  equations  and  thus 
inserted  into  the  equations  of  motion  to  capture  three-dimensional 
induction  effects.  The  equivalent  Theodorsen  function  was  constructed 
with  a  convenient  fractional  calculus  form  to  facilitate  the  capturing 
of  the  frequency-dependent  properties  of  the  generalized  aerodynamic 
forces.  The  coefficient  of  lift  resulting  from  pure  pitching  motion  was 
ccnputed  using  a  method  pres®,  i  which  collapses  to  the  modified 
Doublet  point  method  (27:348-355).  The  coefficient  of  lift  was  then 
used  to  calculate  values  of  the  equivalent  Theodorsen  function  by 
utilizing  two-dimensional  generalized  aerodynamic  force  equations.  The 
values  of  the  equivalent  Theodorsen  function  were  in  turn  modelled  using 
the  fractional  calculus  model.  The  fractional  calculus  model  fully 
captured  the  behavior  of  the  ccnputed  values  of  the  equivalent 
Theodorsen  function  and  coefficient  of  lift  for  rectangular  wings  with 
aspect  ratios  between  two  and  ten.  The  equivalent  Theodorsen  function 


was  generalized  into  the  carpi ex  s -plane  and  the  Wagner  function  and  the 
Kiissner  function,  time -dependent  response  functions  for  indicial  lift 
and  sharp-edge  gusts,  were  determined  and  compared  to  their  two- 
dimensional  counterparts.  The  success  of  capturing  the  properties  of 
the  coefficient  of  lift  for  pure  pitching  motion  for  rectangular 
airfoils  suggests  further  investigation  of  other  generalized  aerodynamic 
forces  caused  by  other  types  of  downwash  for  more  corplex  wing 
geometries.  The  important  calculations  in  aeroelastic  analysis  are  the 
determination  of  the  airloads  and  the  structural  stiffness  and  damping 
matrices.  The  equivalent  Theodorsen  function  allows  the  use  of  the  more 
realistic  three-dimensional  aerodynamic  loads,  rather  than  the  two- 
dimensional  approximate  loads,  in  control  system  designs.  The 
development  of  the  equivalent  Theodorsen  function  for  finite  span 
airfoils  permits  the  general  response  to  be  written  in  a  mathematically 
tractable  form.  It  also  enables  the  active  coupling  between  the 
automatic  control  system  and  aeroelastic  phenomenon  like  flutter  and 
unsteady  structural  loading  which  is  not  available  presently  in  a 
practical  form. 
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A:  Additional  Model  and  Modified  Kernel  Function 
Comparisons  for  Mach  Nvnfcer  0.5 


Model  results  for  the  M=0.5  case  based  on  equation  (26)  are 
further  presented  in  this  appendix.  The  original  discussion  was 
presented  in  Chapter  IV,  compressible  flow  subsection. 

The  model  used  to  produce  the  fit  in  Figure  17  is 


K *  [  0 . 5 ,  Jc,  0 . 1 , 0 . 05]  = 


e.iiw*-  a  ,(■£«»<  (76) 

l+e3(iic)d*  1  +  07  (Ik) *• 


where 


K,  =1.9991 
0,  =0.534 
03=0.116 
@5=1.24 
@7=0.331 


@3=0.523 
04=1.62 
0g=O .  902 
08=1.34 


The  error  is  2. 2x10'*  and  2. 2x10' J  for  the  real  and  imaginary  parts. 
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REDUCED  FREQUENCY 

Figure  17:  Model  of  K*[0.5,1.0,0.05] 
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□  a  a  a  a 


Model  of  K* 

K]  Real  calculated 
K*  Imog  calculated 


REDUCED  FREQUENCY 

Figure  18:  Model  of  K?[0.5,1.0,1.0] 


The  model  for  this  particular  case  in  Figure  18  is  given  by 


K*  [0.5, *,0.1,0.05] 


e.cijc)**  _  e5(iic)*«  (77) 

1  +  03  (lk)$*  l+e7(i*)#* 


where 


K.  =1.7559 
8i  =0.205 
8, =0.302 
8,=1. 87 
87=0.808 


85=0.600 
8i=1.43 
8, =0.954 
8,=1.06 


Errors  are  2. 4x10' J  and  7. 9x10' 3  for  the  real  and  iTOjinary  parts. 
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Additional  sairples  of  the  modified  kernel  function  models  for 
other  regions  in  incompressible  flow  are  presented  in  this  appendix. 

The  original  discussion  is  contained  in  the  inconpressible  subsection  of 
Chapter  IV.  The  form  of  model  was  given  in  equation  (30). 

The  model  used  for  Figure  19  is 

JC*  [ic,  -0.1f0.5]  *  (78) 

1  +  03  {Ik) e* 


where 

K«*=0 . 804 

0i=O.353  9j=0. 93 

0j =0.447  84=0.876 


producing  errors  of  2.3xlO'3  and  1.8x10' 3  respectively. 
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REDUCED  FREQUENCY 

Figure  19:  Model  of  K?[0,-0.1,0.5] 
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ooooo 


Model  of  K* 

K*  Reol  colculoted 
K*  Imog  colculoted 


Appendix  C:  Model  Parameters  for  Incompressible  Case 


Additional  figures  describing  the  variation  of  the  different 
parameters  0j  through  0^  of  the  model  for  the  modified  kernel  function 
in  incompressible  flow  is  further  presented  in  this  appendix.  Other 
figures  were  presented  along  with  the  initial  discussion  in  Chapter  IV. 
The  four-parameter  model  is  given  in  equation  (30). 
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MODEL  PARAMETER 


MODEL  PARAMETER 


CM 


Figure  26: 


CHORDWISE  VARIABLE 

Parameter  2  Dependence  on  Xg 
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MODEL  PARAMETER 


-1.0 


CHORDWISE  VARIABLE 


Figure  27:  Parameter  3  Dependence  on  Xg 
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Appendix  D:  Lift  Coefficient  Model  at  Various  Aspect  Ratios 


The  coefficients  of  lift  for  rectangular  wings  of  various  aspect 
ratios  between  two  and  ten  are  shown  modelled  in  this  appendix.  The 
initial  discussion  and  development  of  these  models  is  presented  in 
Chapter  VI.  The  general  form  of  the  fractional  calculus  model  is 

«  *|iic  +  2r(Jc)|l  +  -4fjj  (80) 

where  the  equivalent  Theodor sen  function  T(k)  is  defined  as  follows. 

r(iJc)  -  T0  +  fl(i*)>*  -  gr(  ik)  +  f(ik)2  (81) 

0  1  +  bUk)* 

This  appendix  will  present  a  series  of  different  figures  showing  the 
superb  agreement  between  the  model  and  the  calculated  lift  coefficients 
for  various  aspect  ratios.  Each  figure  on  the  following  pages  will  be 
acconv>anied  with  an  equation  showing  the  model  of  T(k)  used  to  produce 
the  fit.  The  error  associated  with  each  plot  is  not  explicitly  stated. 
However,  the  error  an  the  real  parts  was  of  the  order  of  10'3  and  the 
error  on  the  imaginary  part  of  the  order  of  10'*  for  all  the  cases 
shown. 


100 


REDUCEO  FREQUENCY 


Figure  29:  Cj ,  for  Wing  of  Aspect  Ratio  2 
C ^  =  n|iic  +  2r(ic)Jl  + 


T{  Ik)  -  T0  +  -  g{ik)  +  f( 

0  1  ♦  b(lk)  ** 


-  Model  of  Ct^ 

ooooo  Real  calculated 

-  Imaginary  Ct,,**  calculated 


2.6 

jx. 


REDUCED  FREQUENCY 


Figure  31:  Cj ,  for  Wing  of  Aspect  Ratio  4 


where 


«  *  jiJc  +  2r(Je)|l  + 


T(ik)  -  r0  ♦  ■  -  ■  -  ff(ik)  +  £{  ik) 

0  1  +  b(lk)* 


Ti  =0.593 
a=0 . 433 
b=2.49 
g=0.042 
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Appendix  E:  MrKHnnal  Model  Parameters  for 
Equivalent  Theodorsen  Function 

The  variation  of  the  parameters  of  the  fractional  calculus  model 
of  the  coefficient  of  lift  are  presented  in  this  appendix.  The 
parameters  and  the  model  are  given  in  equation  (61). 


Figure  36:  Parameter  b  vs.  Aspect  Ratio 
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Pigure  37:  Parameter  Tfl  vs.  Aspect  Ratio 
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ASPECT  RATIO 


Figure  39:  Parameter  £  vs.  Aspect  Ratio 


on 


Based  Upon  References  32  and  16 

The  wing  is  considered  a  nearly  plane  impenetrable  surface  S 
consistent  with  the  concepts  of  linear  theory.  Let  the  wing  lie  nearly 
in  the  x-y  plane  and  let  it  and  the  x-y-s  coordinate  system  to  which  it 
is  referred  be  assumed  to  move  with  uniform  velocity  V  in  the  negative 
x -direction.  Note  the  positive  z -direct ion  is  defined  opposite  to  that 
in  Chapter  2.  At  the  same  time,  let  each  point  of  the  wing  be  assumed 
to  undergo  small  amplitude  harmonic  translations  2^(x,y,t)  at  circular 
frequency  w  and  let  c  represent  the  speed  of  sound  in  the  mediun. 
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The  problem  for  an  oscillating  wing  consists  in  solving  the  wave 
equation  to  certain  boundary  conditions.  The  wave  equation  in 
rectangular  coordinates  is 


(FI) 


The  dependent  variable  2  in  equation  (FI)  is  regarded  as  an 
acceleration  potential.  The  acceleration  is  directly  related  to  a 
perturbation  pressure  field  and  is  related  to  a  velocity  potential  i. 


S 


d« 

Tt 


*  v 


d« 

13? 


(F2) 


The  boundary  problem  for  the  wing  is  ccrrpl eted  by  calculating  the 
downwash  w(x,y,z,t)=di/dz  associated  with  2.  This  downwash  is  assured 
to  be  harmonic  with  regard  to  time  which  implies  that  both  potentials  2 
and  i  are  harmonic  with  respect  to  time  and  as  shown  in  equation  (P3). 


S  (x.y,  z,  t)  «  eimtR(x,y,z) 
•  {x,y,  z,  t)  ■  eimtf(x,y,z) 


(F3) 


With  these  expressions  for  2  and  1,  equation  (F2)  becomes  independent  of 
time  and  reduces  to  an  ordinary  differential  equation  with  one  dependent 
variable 


S  -  <M> 

dx 

This  equation  can  be  integrated  with  respect  to  x  to  produce  equation 
(F5) . 
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2  (X,y,  z)  exp 


it*  A\ 

V  ) 


dX 


(F5) 


The  boundary  problem  for  the  wing  may  now  be  expressed 
mathematically  in  the  following  manner.  Assuming  simple  harmonic 
motion,  equation  (PI)  becomes 


3a8  A  ys  +  d2S 

Bx2  By 2  Bz2 


(F6) 


To  insure  tangential  flow  at  the  wing  surface,  the  potential  must 
satisfy  the  following  downwash  condition. 


w(x'y)  *  J“)z*u'y) 


Here  w  and  Z,  are  anplitudes  of  velocity  and  displacements  respectively 
and  are  assumed  to  be  known  from  the  motion  of  the  wing.  At  z=0,  the 
pressure  must  be  zero  at  all  points  (x,y)  off  the  wing. 

p--p(  K)am0  (F8) 

The  potential  8  is  allowed  to  be  discontinuous  at  all  points  on  the  wing 
and  the  value  of  p  is  determined  by  the  magnitude  of  the  discontinuity 
in  8  at  the  point.  In  the  neighborhood  of  the  trailing  edge,  p  must  go 
to  zero  to  satisfy  the  Kutta  condition.  One  other  condition,  that  9 
vanish  far  ahead  of  the  wing  is  inherently  satisfied  by  the  condition 
given  in  equation  (PS). 
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The  potential  gg  at  a  point  (x,y,z)  due  to  a  harmonically 
pulsating  doublet  located  in  the  x-y  plane  at  a  point  ((,i|,0)  that 
satisfies  equation  (F6)  is 


exp 


Sr 


[Tz 


ift) 

t  +  (x-l)  - 

J?' 

) 

A _ 

cPa 

J 

(F9) 


where 


R'  »  >J{x-\)i  +  Pi(y~n)a  +  pazs 

The  factor  A  is  a  strength  and  dimensionality  factor  allowing  different 
uses  and  interpretations  of  the  potential  gg.  If  gg  is  considered  an 
acceleration  potential  and  substituted  into  equation  (F5),  the 
corresponding  velocity  obtained  may  be  written  in  the  following  manner. 


•o 


e  v 


e 


exdi 


i  to  | 


MX 

cP2 


where 


R  -  )/X 2  +  P2  (y- 11)  2  +  P2r2 


The  downwash  d<g/dz  associated  with  ig  may  be  written  as 


~TF 


.a2  --^3  r*.  oiib-iSF***) 

a*2*  J-  viTCtT77 


(F13) 
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In  equation  (F13) 


*0  -  x-l 

w  -  <i )/V^a 

i  -  V  (y-»i) 5  +  z2 

With  the  use  o£  equations  (F13)  and  (F14)  and  the  concept  of  solving 
linear  boundary- value  problems  by  means  of  superposition  of  elementary 
solutions  to  the  governing  differential  equation,  the  boundary  value 
problem  presented  here  can  be  reduced  to  an  integral  equation. 


w(x,y) 


d2  r*. 

dz*J-~  ^Xa  +  ra 


(F15) 


Here  S  represents  the  surface  of  the  wing  and  L((,q)  represents  an 
unknown  lift  distribution  or  doublet  strength  on  S.  This  integral 
equation  can  be  reduced  further  to  the  form  given  in  Chapter  2  as  shown 
in  Reference  31. 
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